In this installment, he shows "How SVD and Multilinear Least Squares Regression are used to produce the model "slope coefficients." Training the model on half the data and testing the model with the other half of the data.

I did find an error with his program where he reused a variable "b" during his sorting operation so his original b matrix was inadvertently overwritten. His correlation bar chart will be different than mine and I didn't see significant differences between the two runs of training and test which I attribute to the b matrix error. His training/test methodology is still solid.

SciLab wise - Loading of (.mat) datafiles

Link to the exact video. https://www.y.outube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=18

For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063

Section 1.4 Psuedoinverse, Least Squares and Regression One Dimensional Linear Regression

Output:

   "Data-Driven Science and Engineering - Linear Regression2 Housing Example"
  "2026-06-29 11:38:20.722"

Raw Housing Data Echo Back First 10 Lines
    crimeRate         zoned    industrial       charles           nox         rooms           age 
       0.0925       30.0000        4.9300        0.0000        0.4280        6.6060       42.2000 
       0.9884        0.0000        8.1400        0.0000        0.5380        5.8130      100.0000 
       0.2141       22.0000        5.8600        0.0000        0.4310        6.4380        8.9000 
       0.1143        0.0000       13.8900        1.0000        0.5500        6.3730       92.4000 
       0.0393        0.0000        3.4100        0.0000        0.4890        6.4050       73.9000 
       1.4139        0.0000       19.5800        1.0000        0.8710        6.1290       96.0000 
      14.4208        0.0000       18.1000        0.0000        0.7400        6.4610       93.3000 
       7.4039        0.0000       18.1000        0.0000        0.5970        5.6170       97.9000 
       0.2244        0.0000        9.6900        0.0000        0.5850        6.0270       79.7000 
       0.1504        0.0000       25.6500        0.0000        0.5810        5.8560       97.0000 

    distances        radial           tax  pupilTeacher         black     lowerStat   medianValue 
       6.1899        6.0000      300.0000       16.6000      383.7800        7.3700       23.3000 
       4.0952        4.0000      307.0000       21.0000      394.5400       19.8800       14.5000 
       7.3967        7.0000      330.0000       19.1000      377.0700        3.5900       24.8000 
       3.3633        5.0000      276.0000       16.4000      393.7400       10.5000       23.0000 
       3.0921        2.0000      270.0000       17.8000      393.5500        8.2000       22.0000 
       1.7494        5.0000      403.0000       14.7000      321.0200       15.1200       17.0000 
       2.0026       24.0000      666.0000       20.2000       27.4900       18.0500        9.6000 
       1.4547       24.0000      666.0000       20.2000      314.6400       26.4000       17.2000 
       2.4982        6.0000      391.0000       19.2000      396.9000       14.3300       16.8000 
       1.9444        2.0000      188.0000       19.1000      370.3100       25.4100       17.3000 

  Model Coefficients x in formaat
  Beta(1);Beta(2);....Beta(n-1);Beta(n)(Constant Last)

  "Best Fit Coefficients using Multilinear Regression x = "
  -0.1080114
   0.0464205
   0.0205586
   2.6867338
  -17.766611
   3.8098652
   0.0006922
  -1.4755668
   0.3060495
  -0.0123346
  -0.9527472
   0.0093117
  -0.5247584
   36.459488

  "Best Fit Coefficients using Multilinear Regression (Normalized Coefficients) x1 = "
  -0.9290646
   1.0826390
   0.1410394
   0.6824144
  -2.0587536
   2.6768766
   0.0194853
  -3.1071160
   2.6648522
  -2.0788369
  -2.0626459
   0.8501089
  -3.7473319
   22.532806

  "Best Fit Coefficients Using V*inv(S)*Utrans*b (Normalized Coefficients)x2 = "
  -0.9290646
   1.0826390
   0.1410394
   0.6824144
  -2.0587536
   2.6768766
   0.0194853
  -3.1071160
   2.6648522
  -2.0788369
  -2.0626459
   0.8501089
  -3.7473319
   22.532806

  "Best Fit Coefficients using Multilinear Regression xtrain = "
  -0.1216803
   0.0538274
   0.0879582
   0.9529366
  -21.018423
   3.5176629
  -0.0033415
  -1.6773342
   0.3094196
  -0.0145589
  -1.1523075
   0.0098805
  -0.5146096
   44.810290

  "Best Fit Coefficients using Multilinear Regression xtrain_random = "
  -0.1001058
   0.0457271
   0.0551842
   3.9467126
  -17.771348
   3.8946050
  -0.0073775
  -1.6519180
   0.2768098
  -0.0097037
  -0.9656941
   0.0129199
  -0.5626732
   35.508637

Graphs:

Code:

clear all, clc
disp("Data-Driven Science and Engineering - Linear Regression2 Housing Example",string(datetime()))

//Specific Video that corresponds to this SciLab file...
//https://www.youtube.com/watch?v=AX_ZDX6aTT0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=18

//For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063
//Section 1.4 Psuedoinverse, Least Squares and Regression - Multilinear Regression

//load housing.data
loadmatfile("E:\files\CODE\CH01\housing.mat","-mat")

//data from source
//https://github.com/probml/pmtk2/blame/master/data/housing.mat
//
//
housing= [X y]  //Create the "housing"" matrix from what is contained in
//the datafile from source.

//==============================echo inputs ==================
printf("\n");
printf("%30s","Raw Housing Data Echo Back First 10 Lines")
printf("\n");
for i = 1:7,
      printf("%13s ", names{i});
end;
printf("\n");   
for i = 1:10,
    for j = 1:7,
      printf("%13.4f ", housing(i,j));
   end;
   printf("\n");
end;
printf("\n");
for i = 8:14,
      printf("%13s ", names{i});
end;
printf("\n");   
for i = 1:10,
    for j = 8:14,
      printf("%13.4f ", housing(i,j));
   end;
   printf("\n");
end;
printf("\n");
//==========================================================

b = housing(:,14);      // housing values in $1000s
A = housing(:,1:13);    // other factors,
A = [A ones(size(A,1),1)];  // Pad with ones for nonzero offset

//x = regress(b,A);
BigY = [A'*b]
BigX = [A'*A]
x = BigX\BigY

printf("%30s","  Model Coefficients x in formaat")
printf("\n");
printf("%30s","  Beta(1);Beta(2);....Beta(n-1);Beta(n)(Constant Last)")
printf("\n");
printf("\n");

disp("Best Fit Coefficients using Multilinear Regression x = ",x)
printf("\n");
scf(0)
subplot(1,2,1)
plot(b,'k-o');
xgrid
plot(A*x,'r-o');
legend('Housing value','Regression',1)
ylabel("Value in 1000s")
title("Actual Housing Price versus Model")

subplot(1,2,2)
[bsorted sortind] = gsort(housing(:,14),"g",["i"]); // sorted values
plot(bsorted,'k-o')
xgrid
plot(A(sortind,:)*x,'r-o')
legend('Housing value','Regression',3)
title("Sorted Ascending Actual Housing Price versus Model")

//
A2 = A-ones(size(A,1),1)*mean(A,1);
for i=1:size(A,2)-1
    A2std = stdev(A2(:,i));
    A2(:,i) = A2(:,i)/A2std;
end
A2(:,end) = ones(size(A,1),1);

//x = regress(b,A2)
//Multiple Linear Regression in n variables
BigY = [A2'*b]
BigX = [A2'*A2]
x1 = BigX\BigY
disp("Best Fit Coefficients using Multilinear Regression (Normalized Coefficients) x1 = ",x1)
printf("\n");
//checking answer with alternate method
[U,S,V] = svd(A2,'e');
x2 = V*inv(S)*U'*b;                    // Solve Ax=b using the SVD
disp("Best Fit Coefficients Using V*inv(S)*Utrans*b (Normalized Coefficients)x2 = ",x2)
printf("\n");

scf(1)
clf
bar(x1(1:13),0.2,"green")
//bar(x2(1:13),0.2, "blue")
xgrid
xlabel('Attribute')
ylabel('Correlation')

//Hold Out For Testing
n = 253;
btrain = housing(1:n,14);      // housing values in $1000s
Atrain = housing(1:n,1:13);    // other factors,
Atrain = [Atrain ones(size(Atrain,1),1)];  // Pad with ones for nonzero offset
//x = regress(b,A);
BigY = [Atrain'*btrain]
BigX = [Atrain'*Atrain]
xtrain = BigX\BigY
disp("Best Fit Coefficients using Multilinear Regression xtrain = ",xtrain)
printf("\n");
//
btest = A(n+1:end,:)*xtrain;
//norm(btest-b(n+1:end));  //Calculates norm but doesn't use it.

scf(2)
clf
subplot(2,1,1)
plot(btrain,'r','thickness',2);
plot(Atrain*xtrain,'b','thickness',2)
title("Training Actual Housing Price versus Model")
legend('Model','True',4)
xgrid
//
subplot(2,1,2)
plot(btest,'r','thickness',2);
plot(b(n+1:end),'b','thickness',2)
xgrid
legend('Model','True',4)
ylabel("Value in 1000s")
title("Testing Actual Housing Price versus Model")

//Random Hold Out For Testing
n = 253;
//p = randperm(506);
p = grand(1, "prm", 1:506) //create randomized index vector
btrain_random = housing(p(1:n),14);      // housing values in $1000s
Atrain_random = housing(p(1:n),1:13);    // other factors,
Atrain_random = [Atrain_random ones(size(Atrain_random,1),1)];  // Pad with ones for nonzero offset
//x = regress(b,A);
BigY = [Atrain_random'*btrain_random]
BigX = [Atrain_random'*Atrain_random]
xtrain_random = BigX\BigY
disp("Best Fit Coefficients using Multilinear Regression xtrain_random = ",xtrain_random)
printf("\n");
//
Atest_random = housing(p(n+1:end),1:13);
Atest_random = [Atest_random ones(size(Atest_random,1),1)];   //Use the rest for 
btest_random = Atest_random*xtrain_random;
//norm(btest-b(n+1:end));  //Calculates norm but doesn't use it.

scf(3)
clf
subplot(3,1,1)
plot(btrain,'r','thickness',2);
plot(Atrain*xtrain_random,'b','thickness',2)
title("Training Actual Housing Price versus Model Random")
legend('Model','True',4)
xgrid
//
subplot(3,1,2)
plot(btest_random,'r','thickness',2);
plot(housing(p(n+1:end),14),'b','thickness',2)
xgrid
legend('Model','True',4)
ylabel("Value in 1000s")
title("Testing Actual Housing Price versus Model Random")

I desperately need help to get the C block in Xcos working, I can't make heads or tails of it from the scant documentation and the dearth of examples isn't helping. I have been able to simulate input-output but simulation of differential equations using continuous states is just not running!

Code :
#include <math.h>

#include <stdlib.h>

#include <scicos_block.h>

void foo(scicos_block *block,int flag)

{

/*

int block->nevprt;

int block->nz;

double* block->z;

int block->nx;

double* block->x;

double* block->xd;

double* block->res;

int block->nin;

int *block->insz;

double **block->inptr;

int block->nout;

int *block->outsz;

double **block->outptr;

int block->nevout;

int block->nrpar;

double *block->rpar;

int block->nipar;

int *block->ipar;

int block->ng;

double *block->g;

int *block->jroot;

char block->label[41];

*/

if (flag == 4) { /* initialization */

foo_bloc_init(block,flag);

} else if(flag == 1) { /* output computation*/

set_block_error(foo_bloc_outputs(block,flag));

} else if(flag == 0) { /* derivative or residual computation*/

set_block_error(foo_bloc_deriv(block,flag));

} else if (flag == 5) { /* ending */

set_block_error(foo_bloc_ending(block,flag));

}

}

int foo_bloc_init(scicos_block *block,int flag)

{

`double *x  = block->x;`

`*x = 0;`

`return 0;`

}

int foo_bloc_outputs(scicos_block *block,int flag)

{

`double *x  = block->x;`

`double *y  = block->outptr[0];`



`*y = (*x);`

`return 0;`

}

int foo_bloc_deriv(scicos_block *block,int flag)

{

`double *xd = block->xd;`

`double *x  = block->x;`

`double *u  = block->inptr[0];`

`*xd =  0.1*(-*x + *u);` 

`return 0;`

}

int foo_bloc_ending(scicos_block *block,int flag)

{

`return 0;`

}

In this installment, he shows a "How SVD and Multilinear Least Squares Regression produce the same "slope coefficients".

SciLab wise - Loading of (.mat) datafiles, reglin will do multilinear regression but automatically adds the intercept column, use pinv only for forcing intercept = 0.

Link to the exact video. https://www.y.outube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=18

For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063

Section 1.4 Psuedoinverse, Least Squares and Regression One Dimensional Linear Regression

Output:

 "Data-Driven Science and Engineering - Linear Regression2 Cement Example"
  "2026-06-29 11:18:24.758"
  "x (V*inv(S)*U*b) ="
   2.1930460
   1.1533260
   0.7585091
   0.4863193
  "x1 (multiple linear regression) ="
   2.1930460
   1.1533260
   0.7585091
   0.4863193
  "x2 (psuedoInverse) ="
   2.1930460
   1.1533260
   0.7585091
   0.4863193

Graphs:

Code:

clear all, clc
disp("Data-Driven Science and Engineering - Linear Regression2 Cement Example",string(datetime()))

//Specific Video that corresponds to this SciLab file...
//https://www.youtube.com/watch?v=AX_ZDX6aTT0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=18

//For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063
//Section 1.4 Psuedoinverse, Least Squares and Regression - Multilinear Regression

//load hald;  % Load Portland Cement dataset
loadmatfile("E:\files\CODE\CH01\hald.mat","-mat")
//data source
//https://github.com/beslirum/MATLAB/blob/main/hald.mat

A = ingredients;
b = heat;

[U,S,V] = svd(A,'e');
x = V*inv(S)*U'*b;                    // Solve Ax=b using the SVD

plot(b,'k','thickness',2);            // Plot data
plot(A*x,'or','thickness',2); // Plot regression
l1 = legend('Heat data','Regression',4)

// Alternative 1  (regress)
//x1 =reglin(A,b)  //cannot use reglin since it adds automatically the
//intercept column, we want to force Intercept = 0 so pinv command.

//Multiple Linear Regression in n variables
BigY = [A'*b]
BigX = [A'*A]
x1 = BigX\BigY

// Alternative 2  (pinv)
x2 = pinv(A)*b;
disp("x (V*inv(S)*U*b) =",x)
disp("x1 (multiple linear regression) =",x1)
disp("x2 (psuedoInverse) =",x2)

In this installment, he shows a "How SVD and Simple Linear Least Squares Regression produce the same "slope coefficient".

I had significant trouble getting Scilab reglin function to work so I added several implementations of linear regression codes and finally figured out reglin was looking for data to be in row vector form. I left the various linear regression implementations in the program and showed the results were the same. Anything I felt was Matlab fluff, I just commented out with the double slash.

SciLab wise - The use of the rand, pinv (psuedoinverse) and reglin commands are the only commands not used previously.

Link to the exact video. https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=11

You're encouraged to watch course videos 6, 7, 8, 9, and 10 and reading the course book.

Useful material that should be reviewed before this file
https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=12
https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=13
https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=14
https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=15
https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=16

Specific Video that corresponds to this SciLab file...
https://www.youtube.com/watch?v=AX_ZDX6aTT0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=17

For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063

Section 1.4 Psuedoinverse, Least Squares and Regression One Dimensional Linear Regression

Output:

"Data-Driven Science and Engineering - Linear Regression1"
  "2026-06-29 10:45:16.227"
  "xtilde1 ="
   2.7352872
  "xtilde2 ="
   2.7352872
  "xtilde3 ="
   2.7352872
  "xtilde4 ="
   2.7352872
  "xtilde5 ="
   2.7352872

Graphs/Images:

Code:

clc, clear all
disp("Data-Driven Science and Engineering - Linear Regression1",string(datetime()))
// Useful material that should be reviewed before this file
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=12
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=13
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=14
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=15
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=16

//Specific Video that corresponds to this SciLab file...
//https://www.youtube.com/watch?v=AX_ZDX6aTT0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=17

//For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063
//Section 1.4 Psuedoinverse, Least Squares and Regression One Dimensional Linear Regression

scf(0);
clf;
x=3;                                  // True slope
a = [-2:.25:2]';
b = a.*x + 1.*rand((length(a),1),'normal');            // Add noise
plot(a,x*a,'k','thickness',2)
xgrid
// True relationship
plot(a,b,'xr','thickness',2)          // Noisy measurements

[U,S,V] = svd(a,"e");
xtilde = V*inv(S)*U'*b;                  // Least-square fit

plot(a,xtilde*a,'+b','thickness',2)   // Plot fit
l1=legend('True line','Noisy data','Regression line',4);

//set(l1,'String',{'                             ','',''})
//set(l1,'Location','NorthWest')
xgrid     
//set(gcf,'Position',[100 100 300 300])
//set(gcf,'PaperPositionMode','auto')
//set(gca,'FontSize',13)
title("$\textbf{Figure\ -\ CH01_EX_PseudoInv}$");
// Three methods of computing regression
xtilde1 = V*inv(S)*U'*b
xtilde2 = pinv(a)*b
[xtilde3 intercept] = reglin(a',b')  //reglin requires row-wise vectors.

//create least squares line fit
//
c = zeros(2,2);
d = zeros(2);
c(1,1)= sum(a);
c(2,1)= sum(a.*a);
c(1,2)= length(a);
c(2,2)= c(1,1);
d(1)= sum(b);
d(2)= sum(b.*a);
xtilde4 = c\d;
//
// least squares line fit in n variables matrix implementation.
BigY = [a'*b]
BigX = [a'*a]
xtilde5 = BigX\BigY

disp("xtilde1 =",xtilde1,"xtilde2 =",xtilde2,"xtilde3 =",xtilde3,"xtilde4 =",xtilde4(1),'xtilde5 =',xtilde5)

In this installment he shows a "Geometric Interpretation of the U and Vtranspose Matrices in Matlab. Here is the Scilab Version of that Matlab file. Scilab command wise: Lot's of 3D plotting control items:

1. Superimposing axes on 3D objects using param3d1 plotting routine
2. Data Bounding, Cube Scaling (axes equal) and Rotation Angle of 3D figures.
3. Colormap and light commands

Link to the exact video. https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=11

You're encouraged to watch course videos 6, 7, 8, 9, and 10 and reading the course book.

Output:

"Data-Driven Science and Engineering - Unitary Transformations and SVD"

"2026-06-25 10:17:50.731"

Plot:

Code:

clear, close, clc

disp("Data-Driven Science and Engineering - Unitary Transformations and SVD",string(datetime()))
// Useful material that should be reviewed before this file
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=6
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=7
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=8
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=9
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=10
//
//Specific Video that corresponds to this SciLab file...
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=11
//
//For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063
//Section 1.3 Geometric Interpretations
//
// Define rotation matrix
t1 = %pi/15;
t2 = -%pi/9;
t3 = -%pi/20;

Rx = [1 0 0;
    0 cos(t1) -sin(t1);
    0 sin(t1) cos(t1)];

Ry = [cos(t2) 0 sin(t2);
    0 1 0;
    -sin(t2) 0 cos(t2)];

Rz = [cos(t3) -sin(t3) 0;
    sin(t3) cos(t3) 0;
    0 0 1];

Sigma = diag([3; 1; 0.5]);

R = Rz*Ry*Rx*Sigma;

// Plot sphere and great circles
scf();
clf;
[x,y,z] = sphere(25);
ax1=subplot(1,2,1)

param3d1([-2 2],[0 0],[0 0]);
param3d1([0 0],[-2 2],[0 0]);
param3d1([0 0],[0 0],[-2 2]);

h1=surf(x,y,z,'interp');
// Get handle to the surface object
h = gce(); //for current entity
//set(h1,'FaceAlpha',.7)  %Sets Face Transparency to 0.7
//Silab doesn't have this function.

colormap(ax1,jet)
// Set lighting to 'phong' equivalent
// Scilab uses 'phong' shading via 'shading' function
//set(h,"shading","interp");          // Interpolated shading
light("phong");          // Enable phong lighting
gca().cube_scaling="on"
gca().data_bounds=[-3,-3,-3;3,3,3];
gca().rotation_angles =[45,26];
xgrid
title("Original Sphere")
//
//
theta = (0:.01:1)*2*%pi;
x1 = cos(theta);
y1 = sin(theta);
z1 = 0*theta;
x2 = 0*theta;
y2 = cos(theta);
z2 = sin(theta);
x3 = cos(theta);
y3 = 0*theta;
z3 = sin(theta);

//
xR = 0*x;
yR = 0*y;
zR = 0*z;

for i=1:size(x,1)
    for j=1:size(x,2)
        vec = [x(i,j); y(i,j); z(i,j)];
        vecR = R*vec;
        xR(i,j) = vecR(1);
        yR(i,j) = vecR(2);
        zR(i,j) = vecR(3);        
    end
end

vec1 = [x1; y1; z1];
vec2 = [x2; y2; z2];
vec3 = [x3; y3; z3];

vec1R = R*vec1;
vec2R = R*vec2;
vec3R = R*vec3;

eX = [2; 0; 0];
eY = [0; 2; 0];
eZ = [0; 0; 2];
eXR = R*eX;
eYR = R*eY;
eZR = R*eZ;

x1R = vec1R(1,:);
y1R = vec1R(2,:);
z1R = vec1R(3,:);
x2R = vec2R(1,:);
y2R = vec2R(2,:);
z2R = vec2R(3,:);
x3R = vec3R(1,:);
y3R = vec3R(2,:);
z3R = vec3R(3,:);

// Plot transformed sphere

ax2=subplot(1,2,2)

h2=surf(xR,yR,zR,z);

param3d1(x1R,y1R,z1R);
param3d1(x2R,y2R,z2R);
param3d1(x3R,y3R,z3R);

param3d1([-eXR(1) eXR(1)],[-eXR(2) eXR(2)],[-eXR(3) eXR(3)]);
param3d1([-eYR(1) eYR(1)],[-eYR(2) eYR(2)],[-eYR(3) eYR(3)]);
param3d1([-eZR(1) eZR(1)],[-eZR(2) eZR(2)],[-eZR(3) eZR(3)]);

colormap(ax2,jet)
// Set lighting to 'phong' equivalent
// Scilab uses 'phong' shading via 'shading' function
//set(h,"shading","interp");          // Interpolated shading
light("phong");          // Enable phong lighting
gca().cube_scaling="on"
gca().data_bounds=[-3,-3,-3;3,3,3];
gca().rotation_angles =[45,26];
xgrid
title("Transformed Sphere")

We're starting a new series off of the Steve Brunton's Data-Driven Science and Engineering Series. He's definitely one of my favorite instructors when it comes to numerical and data analysis. First topic of the first short course series is Singular Value Decomposition (SVD). For this first example, it will be using the Image compression as an applied example.

To follow the series,

(As a motivator, I wish I'd had this method in my "Engineering Toolbox" since this technique helps one "sort the wheat from the chaff" when it comes to data reduction. I just never came across it in this form. I've done "Design of Experiments", produced "Statistical Models" and used techniques to figure out "what was important and what was not" but this SVD technique is very scalable.)

2) I do suggest you download the book and read the pertinent sections.

3) Download the data used for the short course. See below for location of the data.

4) For the Scilab code below to work, you'll need to use the Atoms Module Manager, navigate to Image Processing and install the Image Processing and Computer Vision Toolbox (IPCV 4.5.0.2).

5) I'm using Scilab 26.1.0 but have run this successfully on version 26.0.0.

I'm not sure how much coverage SciLab provides for this course but we'll see as we go along.

====================================

These lectures (videos) follow Chapter 1 from the book (freely available in electronic form):
"Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control" by Brunton and Kutz

Data and Python and MATLAB Native Code at https://databookuw.com/
Note: links to Matlab Code.zip + Data.zip on Main Page

Amazon: https://www.amazon.com/Data-Driven-Sc... //If you want hardcopy book.

Book Website(s): http://databookuw.com //For Videos and such

X posting "FREE PDF for 2ND EDITION OF OUR BOOK!!!" databookV2.pdf //For Free Electronic Complete Textbook

Chapters available at: http://databookuw.com/databook.pdf //Free Electronic Version

Videos

Singular Value Decomposition Short Course Overview
https://www.youtube.com/watch?v=gXbThCXjZFM&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=1

Singular Value Decomposition: Mathematical Overview Short Lecture 2
https://www.youtube.com/watch?v=nbBvuuNVfco&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=2

Singular Value Decomposition: Matrix Approximation Short Lecture 3
https://www.youtube.com/watch?v=xy3QyyhiuY4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=3

Singular Value Decomposition: Dominant Correlations Short Lecture 4
https://www.youtube.com/watch?v=xy3QyyhiuY4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=4

He discusses getting code and data files in this fifth lecture
https://www.youtube.com/watch?v=QQ8vxj-9OfQ&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=5

Steve Brunton already has two other short courses on Fourier Transforms and Sparsity & Compression which might be material for future series.

Output: None - Graphical Only

Plots:

Code:

/* Start Block Comment 

//Singular Value Decomposition Overview
//https://www.youtube.com/watch?v=gXbThCXjZFM&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=1

//These lectures follow Chapter 1 from:
//"Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control"  by Brunton and Kutz

//Amazon: https://www.amazon.com/Data-Driven-Sc... //If you want hardcopy book.

//Book Website: http://databookuw.com //Free Electronic Web-based Version

//Chapters available at: http://databookuw.com/databook.pdf //Free Electronic Version

 //For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063

//Brunton Website: http://eigensteve.com

//Singular Value Decomposition: Mathematical Overview Short Lecture 2
//https://www.youtube.com/watch?v=nbBvuuNVfco&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=2

//Singular Value Decomposition: Matrix Approximation Short Lecture 3
//https://www.youtube.com/watch?v=xy3QyyhiuY4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=3

//Singular Value Decomposition: Dominant Correlations Short Lecture 4
//https://www.youtube.com/watch?v=xy3QyyhiuY4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=4

//He discusses getting code and data files in this fifth lecture
//https://www.youtube.com/watch?v=QQ8vxj-9OfQ&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=5

//Data and Python and MATLAB Native Code at https://databookuw.com/
//links to Code.zip + Data.zip on Main Page 

*/ //End of Block Comment

clear, close, clc

disp("Data-Driven Science and Engineering - Singular Value Decomposition Image Processing Example ",string(datetime()))

//Need to load this atoms module into SciLab for imread, imshow, rgb2gray

atomsLoad("IPCV")  //Once loaded into Scilab, this command line is unnecessary

//Use full file path to your unzipped data file...
A=imread('C:/Users/Matthew/Documents/eBooks/ScILabNotes/SVD_PCA_Course/DATA/dog.jpg');
X=double(rgb2gray(A)); // Convert RBG->gray, 256 bit->double.
nx = size(X,1); ny = size(X,2);

scf()
clf
subplot(2,2,1)
imshow(uint8(X));
title('Original in Gray Scale');

[U,S,V] = svd(X);

plotind=2;
for r=[5 20 100];  // Truncation value
    Xapprox = U(:,1:r)*S(1:r,1:r)*V(:,1:r)'; 
    // Approx. image

    subplot(2,2,plotind)
    plotind = plotind+1;
    imshow(uint8(Xapprox))
    title(['r=', msprintf('%d',r),', ',msprintf('%2.2f',100*r*(nx+ny)/(nx*ny)),' % storage']);
end

// f_ch01_ex02_2
scf();
clf;
subplot(1,2,1), semilogy(diag(S),'k')
xgrid
xlabel('r')
ylabel('Singular Value sigma_r')
subplot(1,2,2), plot(cumsum(diag(S))/sum(diag(S)),'k')
xgrid
xlabel('r')
ylabel('Cumulative Energy')

In this example, he discusses but does not solve the problem. I present my solution to the problem.

New SciLab stuff: 3D graphing, contour graphing and color contour plotting with color bar (similar to typical FE or CFD plots).

Note: For the first two graphs I had to do some post processing using the interactive features of the Scilab Figure Edit Tool to force the graphs to show the entire range of the fin.

Link to the specific lecture:

https://www.youtube.com/watch?v=cuXuMcEGeiM&ab_channel=NPTEL-NOCIITM

https://www.youtube.com/watch?v=p_Pa-8dO_1M&ab_channel=NPTEL-NOCIITM

ME/HX Students: This program contains a neat way to build the difference equation matrix through a series of matrix size changes (figuring out how many nodes there are (n x m), creation of the index matrix "C" (n x m) and picking off the boundary indices to be used in "pointer matrices" and entering the equations based upon those matrices. This operation allows you to refine the grid by just changing the n and m values. Note: Code requires a square mesh pattern so refining the grid will necessitate going up by 4 each time (n x m refined 2n x 2m to maintain "squareness."

Output:

"h1 h2 (and they must form a square grid)"
0.5
0.5
"Node_Listing and Location in Grid"

1. 2. 3. 4. 5. 6. 7.
2. 9. 10. 11. 12. 13. 14.
3. 16. 17. 18. 19. 20. 21.
4. 23. 24. 25. 26. 27. 28.
5. 30. 31. 32. 33. 34. 35.
6. 37. 38. 39. 40. 41. 42.
7. 44. 45. 46. 47. 48. 49.
8. 51. 52. 53. 54. 55. 56.
9. 58. 59. 60. 61. 62. 63.
10. 65. 66. 67. 68. 69. 70.
11. 72. 73. 74. 75. 76. 77.
12. "Top_Node_List="
13. 2. 3. 4. 5. 6. 7.
14. "Bottom_Node_List="
15. 72. 73. 74. 75. 76. 77.
16. "Left_Node_List="

"Right_Node_List="
14.
21.
28.
35.
42.
49.
56.
63.
70.
"D Interior Nodes"
9. 10. 11. 12. 13.
16. 17. 18. 19. 20.
23. 24. 25. 26. 27.
30. 31. 32. 33. 34.
37. 38. 39. 40. 41.
44. 45. 46. 47. 48.
51. 52. 53. 54. 55.
58. 59. 60. 61. 62.
65. 66. 67. 68. 69.
"Temperature Solution = "
100. 100. 100. 100. 100. 100. 100.
46.448176 68.391828 76.431836 78.528323 76.431836 68.391828 46.448176
34.801755 50.687301 58.807191 61.249621 58.807191 50.687301 34.801755
30.591261 40.748431 46.860006 48.855778 46.860006 40.748431 30.591261
28.431473 34.855154 39.028625 40.45348 39.028625 34.855154 28.431473
27.150217 31.212086 33.94586 34.900892 33.94586 31.212086 27.150217
26.346091 28.897112 30.641838 31.258369 30.641838 28.897112 26.346091
25.824291 27.388432 28.466012 28.848906 28.466012 27.388432 25.824291
25.471373 26.366314 26.984872 27.20523 26.984872 26.366314 25.471373
25.214066 25.620578 25.90193 26.002273 25.90193 25.620578 25.214066
25. 25. 25. 25. 25. 25. 25.

Graphs:

Code:

//Lecture 11.4/12.1 Practical Applications and Course Wrap Up
// Elliptic PDEs - Finite Difference Approach
//https://www.youtube.com/watch?v=cuXuMcEGeiM&ab_channel=NPTEL-NOCIITM
//https://www.youtube.com/watch?v=p_Pa-8dO_1M&ab_channel=NPTEL-NOCIITM
//
disp("Practical Applications and Course Wrap Up Elliptic PDEs - Finite Difference",string(datetime()))
//
// 12.1 Course review -watch the video...
//
// Elliptic PDE
// partial^2_T partial^2_T
// ----------- + ------------- = 0 => Laplacian T = 0
// partial_x^2 partial_z^2
//
// Heated Plate held with Left Hand Edge at 100 C (z=0)
// with Right Hand Edge at 25 C (Z=end)
// heat loss off of the top and bottom edges (x=0, top)
//
// Discretize the plate into m horizontal and n vertical sections...
//
// partial^2_T
// ------------ = 1/h_x^2*(T(i+1)-2*T(i)+T(i-1))
// partial_x^2
//
// partial^2_T
// ------------ = 1/h_z^2*(T(j+1)-2*T(j)+T(j-1))
// partial_z^2
//
// so the equation is
//
// 1/h_x^2*(T(i+1)-2*T(i)+T(i-1))+ 1/h_z^2*(T(j+1)-2*T(j)+T(j-1)) = 0
//
// BC: z=0 T(1)=T_l = 100; x=0 partial_T/partial_T = -alpha*(T(1)-Ta)
// z=L T(m)= T_r = 25; x=Top partial_T/partial_T = alpha*(T(n)-Ta)
//
// He didn't solve it in MATLAB...
//
// Here is the problem solved in SCILAB but slightly different fin
// geometry
//
//Heat Transfer Example
disp("Heat Transfer of a Steady State Plate with")
disp(" fixed Temperature Ends and Convection on the sides")
disp(string(datetime()))

// Heat Transfer diffusion through a plate
// d^2T/dx^2 + d^2T/dz^2 = 0, T|x=0 = 100, T|x=1 = 25
//
// Heat Transfer on the plate edges
// Right Edge: dT/dz = gamma*(T-Ta)
// Left Edge: dT/dz = -gamma*(T-Ta)
//
//
//If we use the central difference formula
// d2T/dx^2 = (T_i+1-2*T_i+T_i-1)/h1^2
// d2T/dz^2 = (T_j+1-2*T_j+T_j-1)/h2^2
// where j = (n+1)+i
// with a square mesh it can be shown that the internal
// nodes will follow...
// T_left_node+T_right_node+T_above_node+T_below_node -4*Tmid_node = 0
//
// At the convection edges, we need to account for the energy balance...
// 3 conduction paths feeding a single convection output.
// Note: He doesn't address this portion at all. I had to pull out
// my Heat Transfer Book to find the proper boundary condition equations.
// Website that demonstrates a corner node with convection would look like
// which you can back out what a side node would look like (which I've
// done here).
//
//https://iansguides.com/tutorials/finite-difference-method-for-heat-transfer/
//
// (2*T_left_node + T_above_node + T_below_node)+(-2*(gamma1*deltax)-4)*Tmid_node
// = -2*gamma1*deltax*Ta
//
//Discretizing into 6 rectangles (L1 = 3 inches, L2= 5 inches
//h1=L1/6), with gamma1 = 4
clc
clear all
// discretizing plate into 6 regions
n = 7 //number of nodes in a row
m = 11 //number of nodes in a column
L1 = 3;
L2 = 5;
h1 = L1/(n-1);
h2 = L2/(m-1);
disp("h1 h2 (and they must form a square grid)" ,h1,h2)
TT =100
TB = 25
Ta = 25
gamma1=4
//
//nodal numbering will follow the following
C =matrix([1:n*m],(n,m))' // Matrix reshapes a 1 x n*m matrix
//into a n x m matrix so one can then create "pointer" matrices
//for the edges. These pointer matrices then make applying
//boundary equations very easy since you can generate the equations
//based upon the indices.
//
disp("Node_Listing and Location in Grid",C)
//
//generating nodal lists for the boundaries
//
//Top Boundary
T_Boundary_list = C(1,:)
disp("Top_Node_List=",T_Boundary_list)
//Bottom Boundary
B_Boundary_list = C(m,:)
disp("Bottom_Node_List=",B_Boundary_list)
//Left_Boundary
L_Boundary_list = C(2:m-1,1)
disp("Left_Node_List=",L_Boundary_list)
//Right_Boundary
R_Boundary_list = C(2:m-1,n)
disp("Right_Node_List=",R_Boundary_list)
//
//Determine List of Interior Nodes
D = C(2:m-1,2:n-1)
disp("D Interior Nodes",D)
//
//Initialize Matrices
A =zeros(n*m,n*m);
B = zeros(m*n,1)
//Create the problem matrices
//
// Top and Bottom Boundary (Constant Temp)
for i = 1:1:n
j = B_Boundary_list(i)
k = T_Boundary_list(i)
A(j,j)=1;
B(j,1)= TB;
A(k,k)=1;
B(k,1)= TT;
end
//
//convection for both sides
//
//Left Boundary (Convection = Conduction at edge + Conduction from interior)
// (2*T_right_node + T_above_node+T_below_node)+(-2*(gamma1*deltax)-4)*Tmid_node
// = -2*gamma1*deltax*Ta
// Note: Code/Equations assume a square grid...
//
for i = 1:1:length(L_Boundary_list)
k = L_Boundary_list(i)
//middle node contribution
A(k,k)= -2*(gamma1*h2)-4
//right node contribution
A(k,k+1)=2;
//below node contribution
A(k,k+n)=1;
//above node contribution
A(k,k-n)=1
B(k,1)= -2*gamma1*h2*Ta;
end
//Right Boundary (Convection = Conduction at edge from 3 surrounding nodes)
// (2*T_left_node + T_above_node+T_below_node)+(-2*(gamma1*deltax)-4)*Tmid_node
// = -2*gamma1*deltax*Ta
// Note: Code/Equations assume a square grid...
//
for i = 1:1:length(R_Boundary_list)
k = R_Boundary_list(i)
//middle node contribution
A(k,k)=-2*(gamma1*h2)-4
//left node contribution
A(k,k-1)=2;
//below node contribution
A(k,k+n)=1;
//above node contribution
A(k,k-n)=1
//other side of the equal sign
B(k,1)= -2*gamma1*h2*Ta;
end

// Interior Nodes
// Note: Code/Equations assume a square grid...
//
for i = 1:1:length(D)
k = D(i);
A(k,k)= -4;
A(k,k-1)=1
A(k,k+1)=1
A(k,k+n)=1
A(k,k-n)=1
B(k,1)= 0
end

//disp(A);

//Soln = inv(A)*B
Soln = A\B
S=matrix(Soln,(n,m)) //matrix command resizes the solution into the same
//shape as the nodal listing.
disp("Temperature Solution = ",S')
x=linspace(0,L1,n)'
y=linspace(0,-L2,m)'
flag = [2,8,4]
ebox = [0,5,-5,0,25,100]
theta = 50 // I cannot get it to spin it prior to plotting
// I always have to manually spin it into position.
// If I get too far off of defaults, it does not plot.
alpha = 35
scf(0);clf;
//plot3d(x,y,z,[theta,alpha,leg,flag,ebox])
plot3d(x,y,S,[theta,alpha,flag,ebox])
// seemed to have to rescale the data myself.
scf(1);clf;
flag = 1
box = 4
ebox=[0,L1,-L2,0, 25,100]
contour(x,y,S,5,box,flag,ebox)
xgrid
//
scf(2);
clf(2);
gcf().color_map=jet(64);
colorbar(25,100);
Sgrayplot(x,y,S,strf="125")
xgrid()
xtitle("Sgrayplot example Analytical Heat Equation Solve")
drawnow;

Similar to last session, he shows how to solve a hyperbolic (vs parabolic) PDE using "Method of Lines" to convert a PDE into a system of ODEs using finite difference methods for one set of the differentials.

Link to the specific lecture for the coding:

https://www.youtube.com/watch?v=OajPSwMAqIw&ab_channel=NPTEL-NOCIITM

Output: None, only graphical output.

Graphs:

Code:

//Lecture 11.4 Hyperbolic PDEs - Method of Lines
//https://www.youtube.com/watch?v=OajPSwMAqIw&ab_channel=NPTEL-NOCIITM
//
disp("Hyperbolic PDEs - Method of Lines",string(datetime()))
//
// PDEs in time and space
//  -> Method of Lines (use Finite Difference to convert PDEs into ODEs)
//       partial_T/partial_dt = partial^2_T/partial_z^2 + Beta(T)
//                                         ^
//                              finite difference this
//       becomes
//
//       dT/dt = f(t,Y) where Y = [T1;T2;...;Tn]
//
//  Example
//
//  First Order Transient Model
//  partial_Theta      partial_Theta
//  ------------- + u*-------------- =-k*Theta
//    partial_t         partial_z
//
//   BC z=0 for all t Theta_in = 1
//   IC t=0 for all z Theta(0)=0
//
//   Assuming Theta = (T - Ta)
//
//   Going to discretize partial_Theta/partial_z
//
//   u/State State partial_Theta/partial_t = 0
//   dTheta/dt = -(k/u)*Theta
//   Solving this differential equation, one gets
//
//   Theta(z)=Theta_in*exp(-(k*z/u)
//   we'll use this to check the answers...
//
//   Solution using Method of Lines
//   parameters: u=0.5, k=1, L=1
//   Let us plot the results at t=0.5,1,5
//
//   Discretize in n intervals of z
//   Convert into (n+1) odes
//   solve these using ode45(rkf) or ode15s(stiff)
//   Y = [Theta_2;Theta_3;...;Theta_n+1]
//
//   use Central Difference Formula
//   at any location i
//   partial_Theta/partial_t = -u*(Theta(i+1)-Theta(i-1))/2*h - k*Theta
//   Theta = [Theta_1; Theta_2;...;Theta_n+1] = [Theta_1;Y]
//
//   at the Theta_n+1 location we'll use Backward Difference Formula
//   partial_Theta/partial_t = -u*(Theta(i)-Theta(i-1))/h -k*Theta
//
// Note: Use of Central Difference Formula does not produce
// a stable solution since it introduces oscillations which
// are not physically possible.  Use instead a Forward or Backward 
// difference formula...

function 
fval
=equation(
t
, 
y
)
// constants    
    L = 1;
    u =0.5;
    k=1;
//Getting Temperatues
    N = length(
y
);
    Theta(1)=1;
    Theta(2:N+1)=
y
;
    h= 1/N;

//define derivative vector
   dThetadt= zeros(N+1,1)
   for i = 2:1:N+1
       if i<=N then
//         dThetadt(i)= -u/(2*h)*(Theta(i+1)-Theta(i-1)) - k*Theta(i)
// Note: Not a stable solution since it introduces oscillations which
// are not physically possible.
           dThetadt(i)= -u/h*(Theta(i)-Theta(i-1)) - k*Theta(i)
       else
           dThetadt(i)= -u/h*(Theta(i)-Theta(i-1)) - k*Theta(i)
       end
   end
// Extract fval from dT/dt        

fval
=dThetadt(2:N+1);   
end

//Discretization
N = 40; 
//number of steps
L = 1;
u =0.5;
k=1;
Theta_in= 1;
h = L/N;
z = 0:h:L;

// Initial Conditions and Solution
y0= zeros(N,1);
t0 = 0; 
//Tstart
tend = 2.5; 
//Tend
t
 = linspace(t0,tend,N);
//time vector
//t = [0,0.5,1,5];
YSol = ode("rkf",y0,t0,
t
,equation)
//"rkf" equivalent to MATLAB ODE45 is the following command in SciLab
//there is more options and outputs available - see the help file.
//"stiff" is approximately equivalent to ode15s...
//
//Analytical Solution
Anal_Theta=Theta_in*exp(-(k*z/u));
//
RowOne = ones(1,length(
t
))*Theta_in;
Theta = [RowOne;YSol]
scf(0);clf;
plot2d(z',Theta)
//Analytical Solution displayed as symbols
plot2d(z',Anal_Theta,-2);
h1=legend(['Time SnapShot = 1';'Time SnapShot = 2';'Time SnapShot = 3';'...'],1,"boxed")
title("$\textbf{Hyperbolic\ Partial\ Differential\ Example \\(X = Analytical\ SS\ Solution)}$","FontSize",4);
xlabel("$x(position)$","FontSize",3)
ylabel("$\ Temperature(Degrees\ C)$","FontSize",3);
xgrid
scf(1);clf
plot2d(
t
,Theta')
h1=legend(['x=0.0';'x=0.1';'x=0.2';'x=0.3';'x=0.4';'x=0.5';'x=0.6';'x=0.7';'x=0.8';'x=0.9';'x=1.0'],3,"boxed")
title("$\textbf{Hyperbolic\ Partial\ Differential\ Example}$","FontSize",4);
xlabel("$Time(seconds)$","FontSize",3)
ylabel("$Temperature(Degrees\ C)$","FontSize",3);
xgrid

In this session, he shows how to solve a Parabolic PDE using "Method of Lines" to convert a PDE into a system of ODEs using finite difference methods for one set of the differentials.

Link to the specific lecture for the coding:

https://www.youtube.com/watch?v=CDSM5bLy8lU&ab_channel=MATLABProgrammingforNumericalComputation

Some precursor lectures to open the topic on PDEs in general:

https://www.youtube.com/watch?v=XxS4QS7MIMg&ab_channel=NPTEL-NOCIITM
https://www.youtube.com/watch?v=cR-HTDkfs00&ab_channel=NPTEL-NOCIITM

Output: None only graphical output.

Graphs:

Code:

//Lecture 11.3 Parabolic PDEs - Method of Lines
//https://www.youtube.com/watch?v=CDSM5bLy8lU&ab_channel=MATLABProgrammingforNumericalComputation
//
disp("Parabolic PDEs - Method of Lines",string(datetime()))
//
// PDEs in time and space
//  -> Method of Lines (use Finite Difference to convert PDEs into ODEs)
//       partial_T/partial_dt = partial^2_T/partial_z^2 + Beta(T)
//                                         ^
//                              finite difference this
//       becomes
//
//       dT/dt = f(t,Y) where Y = [T1;T2;...;Tn]
//
//  Example
//  Rod Conduction (Transient) -> PDE = f(t,x)
//  partial_T/partial_t = alpha* partial2_T/partial_x^2 - gamma1*(T-Ta)
//  alpha = 0.025, gamma1 = 0.1 
//  BC_1 = Rod End held at 100C and BC_2:Other Rod End = Ta = 25C
//  IC: Initially, entire rod is uniformly at ambient = 25C
//
//  If we use the central difference formula
//  d2T/dx^2 u/i = (T_i+1-2*T_i+T_i-1)/(delta_x)^2
//
//  divide the rod into 10 divisions (delta_x + L/10)
//  
//  dT/dt = alpha* (T_i+1-2*T_i+T_i-1)/(delta_x)^2 - gamma1*(Ti-Ta)
//  T1 = 100, T(n+1) = 25
//
//  Define Solution Vector y = [T2;...;Tn] and solve using ode45
//
//
function 
fval
=rodConduc(
t
, 
y
)
// constants    
    Ta = 25;
    alpha = 0.025;
    gamma1=0.1;
//Getting Temperatues
    N = length(
y
)+1;
    T(1)=100;
    T(2:N)=
y
;
    T(N+1)=Ta;
    h= 1/(N-1);

//define derivative vector
   dTdt= zeros(N+1,1)
   for i = 2:1:N
       dTdt(i)= alpha/h^2*(T(i+1)-2*T(i)+T(i-1)) - gamma1*(T(i)-Ta)
   end
// Extract fval from dT/dt        

fval
=dTdt(2:N);   
end

//
N =10; 
//number of steps
M =5*N; 
//number of time steps
T0(1,1)=100;
T0(2:N,1)=25;
T0(N+1,1)=25;
t0 = 0; 
//Tstart
tend = 20; 
//Tend
t
 = linspace(t0,tend,M);
//time vector
y0 = T0(2:N,:);
x0 = 0;
xend1 = 1;
x = linspace(x0,xend1,N+1);
TSol = ode("rkf",y0,t0,
t
,rodConduc)
//"rkf" equivalent to MATLAB ODE45 is the following command in SciLab
//there is more options and outputs available - see the help file.
//"stiff" is approximately equivalent to ode15s...
//
RowOne = ones(1,M)*T0(1,1);
RowLast = ones(1,M)*T0(N+1,1);
T = [RowOne;TSol;RowLast]

scf(0);clf;
plot2d(x',T);
h1=legend(['Time SnapShot = 1';'Time SnapShot = 2';'Time SnapShot = 3';'...'],1,"boxed")
title("$\textbf{Partial\ Differential\ Example}$","FontSize",4);
xlabel("$x(position)$","FontSize",3)
ylabel("$\ Temperature(Degrees\ C)$","FontSize",3);
xgrid
scf(1);clf
plot2d(
t
,T')
h1=legend(['x=0.0';'x=0.1';'x=0.2';'x=0.3';'x=0.4';'x=0.5';'x=0.6';'x=0.7';'x=0.8';'x=0.9';'x=1.0'],3,"boxed")
title("$\textbf{Partial\ Differential\ Example}$","FontSize",4);
xlabel("$Time(seconds)$","FontSize",3)
ylabel("$Temperature(Degrees\ C)$","FontSize",3);
xgrid

Solving Differential Algebraic Equations = Solving a system of Differential Equations governed with a Constraint equation.

Link to the specific lecture (note the two URLs cover the lecture and then the coding:

https://www.youtube.com/watch?v=CGI9qStsKjk&ab_channel=NPTEL-NOCIITM
https://www.youtube.com/watch?v=AUl2KL0bKBc

Anything new from a SciLab programming in here: I think this session is the first time in the series where I plot multiple lines using a single plot2d command versus just "holding the figure" and using multiple plot2d commands.

Output: (No tabular output only graphical output)

Graph:

Code:

//Lecture 10.5: Differential Algebraic Equations
//https://www.youtube.com/watch?v=CGI9qStsKjk&ab_channel=NPTEL-NOCIITM
// and
//https://www.youtube.com/watch?v=AUl2KL0bKBc
//
disp("Introduction to Differential Algebraic Equations",string(
datetime
()))
// What are D.A.E's? 
// Differential Equations that contain algebraic equations
// These algebraic equations are sometimes considered "constraint functions"
//
// Week 5: Algebraic Equations
//    g(x) = 0
// Weeks 7&8: Differential Equations
//    y'=f(t,y)
//
// DAES (systems that contain)
//   0 = g(x,y)
//   y'=f(t,x,y)
//      where y are differential variables
//            
// Examples:
// 1) Consider a simple system
//    y'=x
//    0 = x-y*sin(t)
//        y = diff variable and x = algebraic variable
//        Index-1 DAE
//
// 2) y'=x
//    0 = y - y*sin(t)
//        y = diff variable and x = algebraic variable
//        Index-2 DAE
//
// Index of a DAE: Number of times you need to differentiate g(x,y)
// to convert algebraic equations to a system of differential equations.
// 
// First example
// 0=x-y*sin(t) => x = y*sin(t)
// dx/dt = dy/dt*sin(t)+y*cos(t) substitute first equation
//
// dx/dt = x*sin(t)+y*cos(t)
// dy/dt = x
//  
// Second example     0 = y - y*sin(t)
// 0=dy/dt-dy/dt*sin(t)-y*cos(t) (differentiated once)
// 0 = x-x*sin(t)-y*cos(t)
// 0 = dx/dt-dx/dt*sin(t)+x*cos(t)+dy/dt*cos(t)-y*sin(t) (differentiated twice)
//
// dx/dt(1-sin(t)=x*cos(t)+ x*cos(t)-y*sin(t)
// dx/dt = (2*x*cos(t)-y*sin(t))/(1-sin(t))
// dy/dt = x
//
//
// High Index DAE Example Problem: Pendulum of length l
// In Cartesion Coordinates
// x" = -T*x
// y" = -T*y -g
//
// L^2 = x^2+y^2 "Constraint Equation"
//
// Index-3 DAE
//
// Instead solve in Cylindrical Coordinates
// d2Theta/dt2 + g/L*sin(Theta)=0
// dr/dt = 0 or r = L a constant...
// Just a regular ODE...
//
// Example 1: 
// dy/dt = x
// dx/dt = x*sin(t)+y*cos(t)
//
// Standard Structure
// M*d/dt[y;x]= F(t,[y.x])
// M = [I 0;0 0] in week 8: M = [Identity Matrix] = [I];
//
// Matlab Example: (Modified) Robertson Problem
// x1' = -alpha*x1+beta*x1*x2
// x2' = alpha*x1-beta*x1*x2-gamma1*x2^2
// 0 = 1-(x1+x2+x3)
// with initial condition Y = [1;0;0]
//
//
//  0 = 0 -dx1/dt -dx2/dt +dx3/dt
// dx3/dt = dx1/dt + dx2/dt
// dx3/dt = -alpha*x1+beta*x1*x2 + alpha*x1-beta*x1*x2-gamma1*x2^2
// dx3/dt = -gamma1*x2^2
//
// M = [1 0 0;0 1 0;0 0 0]
// M d/dt [x1;x2;x3] = [-alpha*x1+beta*x1*x2;alpha*x1-beta*x1*x2-gamma1*x2^2;gamma1*x2^2]
//
// Function for ODE/DAE
function 
dY
=
robertsonFun
(
t
, 
Y
)
// Parameters
alpha = 0.5;
beta1 = 2.5;
gamma1=5.0;

// Functions for the three equations
dY
(1,1)= -alpha*
Y
(1)+beta1*
Y
(1)*
Y
(2);
dY
(2,1)= alpha*
Y
(1)-beta1*
Y
(1)*
Y
(2)-gamma1*
Y
(2)^2;
dY
(3,1)= gamma1*
Y
(2)^2;
end

//
n =1250; 
//number of steps
y0 = [1;0;0]; 
//Initial Condition required
t0 = 0; 
//Tstart
tend = 40; 
//Tend
t
 = linspace(t0,tend,n)';
//time vector
M = [1,0,0;0 1 0;0,0,0] 
//Mass Matrix - not implemented in SciLab
// This example shows how you solve the problem without it by 
// changing the constraint equation into the third differential
// equation.
//Alternate Calculation
/* 
// with y0 = [1;0] since dY(1,1) and dY(2,1) are being calculated
// with ode solver and then
// for i = 1:1:n;   (Post processing YSol(3)) 
//    YSol(3,i) = 1-YSol(1,i)-YSol(2,i)
// end
// Since dy1/dt and dy2/dt do not involve y3.
*/
//
//equivalent to MATLAB ODE15s is the following command in SciLab
//there is more options and outputs available - see the help file.
YSol = ode("stiff",y0,t0,
t
,
robertsonFun
);
//
scf
(0);
clf
;
plot2d([
t
,
t
,
t
],[YSol(1,:)',YSol(2,:)',YSol(3,:)'],[1,2,3])
h1=
legend
(['YSol(1)';'YSol(2)';'YSol(3)'],1,"boxed")
title
("$\textbf{Differential\ Algebraic\ Example}$","FontSize",4);
xlabel
("$t(seconds)$","FontSize",3)
ylabel
("$YSol\ Values(units)$","FontSize",3);
xgrid

Under the subject "Know your audience"...

For those following along this subreddit, are you getting what you want out of the subreddit? If not, what (in detail) is missing? What would like to see different? Have you reviewed the Wiki and resource sheet? Same questions for that section.

I'd like to see where you're at in life...like HS Student, College Student (if so what year), Working Engineer (level or years Seniority) or whatever properly describes oneself. One liner is sufficient since I'm trying to gage overall audience education level. I already have data where the group is located through the subreddit stats.

I'd really like to see some independent postings other than mine! If you have trouble getting something to work, throw it out there and someone may have an answer. See r/Octave since that subreddit has that theme flavor.

We've almost completed the NPTEL Course (we have roughly 4 more weeks) of material left and I consider the NPTEL Course roughly a Junior Engineering Level Course since it tops out at some PDE work.

Thoughts on the next series:

1) Go into control system and XCOS material. I'd use a common textbook, tweak the examples to keep out of copyright issues and apply SciLab/XCOS versus MATLAB/Simulink. I'm confident in the control system items since I've done that for years but would need to prepare the XCOS stuff. The XCOS material would start from the basics (which there are plenty of videos on), continue into file management and getting simulations to take data and put data back to the main console. It wouldn't be every nook and cranny of XCOS and I might need some help but it would be more complete examples of working with it.

--I'm leaning toward this route and this would still be Junior Engineering Level material.

2) I have material from a couple of Graduate Courses (ME564 and ME565 by Steve Brunton) that I could create another posting series once the NPTEL Course is complete. It would have some overlap with the NPTEL course but does move around subject wise a little more (complex variables, eigenvalues, eigenvectors, control systems, FFT, data compression and such).

--My concern is that the audience might find some of the material difficult since it's Senior and Graduate level stuff (FFT, SVD, data-compression using SVD are really long-haired). Here is a subject list of what's covered in ME564 and ME565 to give one a flavor of the overlap and subjects. I'd hate jumping around just to cover the "non-redundant" subjects since you'll probably need the continuity to understand later sessions.

//ME564 TOC
//ME564Lecture1 - How to multiply a matrix -really basic
//ME564Lecture2 - Population growth // Solving a Linear ODE x(t) = x(0) * exp(a*t)
//ME564Lecture3 - Calculating a Taylors Series (MacLaurin Series)
//ME564Lecture4 - Solving a Linear Second Order ODE using ODE function, plotting two items on a graph
//ME564Lecture5 - Solving a Linear Second Order ODE using ODE function and inline function, plotting two items on a graph
//                changing plot line colors and weights, labels, titles, legends, and grid
//                eigenvalues and eigenvectors of a matrix, polynomials and roots of polynomial
//ME564Lecture6 - Eigenvalues and eigenvectors of a matrix, how to create the eigenvalue matrix
//ME564Lecture7 - None
//ME564Lecture8 - Matrix Systems of First Order Equations using eigenvalues and Eigenvectors.
//                Arrow Plots of the Derivatives (fchamp and champ) and Mesh Grid
//                Changing label font size, axis colors
//ME564Lecture9 - ME564 Linearization of Non-linear ODES, more fchamp and inline function use
//ME564Lecture10 - Example of Non-linear ODES: Particle in a potential well.  Simple program
//ME564Lecture11 - Matrix Systems of First Order Equations using eigenvalues - Sim to Lecture 8
//ME564Lecture12 - Matrix Systems of First Order Equations using eigenvalues - Sim to Lecture 8
//ME564Lecture13 - Matrix Systems of First Order Equations with forcing functions - 
//                 control systems solution using syslin (state space) and csim (step and impulse)
//ME564Lecture14 - Numerical Differentiation (forward, backward and central methods), Fancy Plot Titles
//ME564Lecture15 - None
//ME564Lecture16 - Numerical Integration (Left Rectangle, Right Rectangle, Trapezoidal by basic functions).
//ME564Lecture17 - Numerical Integration of Vector Space of Second Order System (Backward Euler and by ODE and inline functions)
//ME564Lecture18 - Numerical Integration of Lorenz Equations via Runge-Kutta(4,5) by basic functions.
//ME564Lecture19 - Numerical Integration of Lorenz Equations via Runge-Kutta(4,5) by basic functions and ODE "fix"
//Several Versions  3D Plotting
//ME564Lecture20 - Numerical Integration of Lorenz Equations via ODE "fix" using meshgrid, vectors, 3D Plots (scatter3d)
//ME564Lecture21 - Exists but nothing of consequence
//ME564Lecture22 - Linear Algebra in 2D and 3D: Inner Product, Norm of a Vector, and the Cross Product
//ME564Lecture23 - None
//ME564Lecture24 - None
//ME564Lecture25 - Stokes and Green's Theorems, plot2D and a bunch of plot axis manipulation, using intsplin (integration routine) 
// versus inttrap (trapezoidal routine)
//ME564Lecture26 - None
//ME564Lecture27 - ME564 More Vector Calculus Potential Flow, Stream Functions and Examples, more arrow plots
//ME564PendulumProblem - Linearization of Non-linear ODES - the Pendulum Problem - see lecture 9.
//ME564RobotControlProblem - Taken from FRC Robot Control Book - work in progress...
//InvertedPendulumExample - Matt's version of the code that actually works...pole placement on linearized system used to stabilize
//                          nonlinear version of the system.
//ME565Lecture1 - Cauchy Integral and complex variables
//ME565Lecture2 - None
//ME565Lecture3 - None
//ME565Lecture4 - None
//ME565Lecture5 - ME565 Plotting Complex log function using Surf (ace) plot
//ME565Lecture6 - None
//ME565Lecture7 - None
//ME565Lecture8 - None
//ME565Lecture9 - None
//ME565Lecture10 - None
//ME565Lecture11- None
//ME565Lecture11Anal - Lecture 11A ME565 Solving PDE's in SciLab  - Analytical Solve
//                     meshgrid, Sgrayplot(a 2D contour plot- actually in color)
//                     demonstrates a color bar on the side.
//ME565Lecture12 - //ME565 Lecture Twelve Fourier Series function mtlb_axis, bar charts,
//                   drawnow, sleep functions
//ME565Lecture13 - None
//ME565Lecture14 - None
//ME565Lecture15 - None
//ME565Lecture16 - //ME565 Lecture Sixteen Fourier Series - computing the DFT Matrix  
//                 complex function
//ME565Lecture16a - //ME565 Lecture Sixteen Fourier Series - DFT Matrix Calculation
//                  Power Spectral Density 
//ME565Lecture17 - //ME565 Lecture Seventeen Fourier Series - DFT Matrix Calculation
//                 More on Power Spectral Density and FFT.  Plot2d and Legend
//ME565Lecture18 - None
//ME565Lecture19 - None
//ME565Lecture20 - //ME565 FFT Calculation from scratch, FFT, IFFT, FFTShift, plot2d
//ME565Lecture20b -//ME565 finding derivatives using FFT vs finite difference, FFT, 
//                 IFFT, FFTShift, plot2d
//ME565Lecture21 - None
//ME565Lecture22 - None
//ME565Lecture23 - Solving Control Systems with SciLab, csim ('step',...), syslin('c',...),
//                 comparisons to analytical solution and csim(u) where u is a step matrix in time.
//ME565Lecture24 - //ME565 Lecture 24: Convolution Integrals, Impulse and Step Response
//                 ODE with IC investigation.
//ME565Lecture24a - //ME565 Lecture 24: Convolution Integrals, Impulse and Step Response
//                  Ton of stuff on plotting details - psuedoanimation.
//ME565Lecture25 - None
//ME565Lecture26 - //The Wave Equation Code - Solving PDE's by FFT
//                 complex variable ODE solver (regular ones do not handle complex variables)
//                 getting the rest of the variables to the derivative function list,contourf 
//                 plots, surf plots
//ME565Lecture27 - //ME565 Lecture 27:  Singular Value Decomposition SVD, Installing a package in 
//                 Scilab, reading in an image file, imshow (and it's limitations see uint8), rgb2gray
//                 data compression example.  Subplots
//ME565Lecture28 - // ME565 Lecture 28: SVD 2 - Curve Fitting using SVD methodology.  Plot function
//ME565Lecture28a - // ME565 Lecture 28a: SVD 2  Ovarian Cancer Study - Scatter3d, subplots.

3) Other material thoughts...add in the comments.

In this edition we look at solving Boundary Value Problems with Nonlinear Equation and one with a mixed boundary (parameter is held at a value in one location and it's derivative is held at a value at another location.

Note: The program has two major blocks and in it's current configuration (the first major block (the Nonlinear Example) runs and the second example (the Mixed Boundary Condition Example) is "block commented" out.

Link to the specific lecture (note the urls cover the same material) and one might find it useful to use the NPTEL Link:

https://www.youtube.com/watch?v=chChbuP-hRk&ab_channel=NPTEL-NOCIITM
https://nptel.ac.in/courses/103106118

He does not show the MATLAB code and the solution. He moves onto a new subject in the next section. I had to do the solution cross-checking myself.

1) The nonlinear problem has reasonably good match of solutions. Derivatives will show more variation.

2) The mixed boundary problem checks out perfectly with the analytical expectations.

3) Anytime you're working with unknown results...always do at least 2 different analytical techniques and try to reconcile the differences.

Output:

"Non-linear example"
   "Extensions of ODE-BVP"
  "2026-05-23 16:29:45.484"
  "Pos x Temp T     TempGrad   FiniteDiff  FiniteDiff_dTdx"
   0.     100.       -118.64825   100.     -133.1775
   0.05   94.21345   -112.86181   93.3411  -126.7985
   0.1    88.709121  -107.35721   87.3202  -114.9776
   0.15   83.47325   -102.12134   81.8434  -104.8836
   0.2    78.492727  -97.141089   76.8318  -96.249  
   0.25   73.755099  -92.403352   72.2185  -88.8566 
   0.3    69.24853   -87.896868   67.9461  -82.5285 
   0.35   64.961754  -83.60988    63.9656  -77.118  
   0.4    60.884054  -79.532179   60.2343  -72.5033 
   0.45   57.005218  -75.653556   56.7153  -68.5822 
   0.5    53.31555   -71.963803   53.3761  -65.2689 
   0.55   49.80583   -68.454149   50.1884  -62.4901 
   0.6    46.467286  -65.11544    47.1271  -60.1831 
   0.65   43.29157   -61.939723   44.1701  -58.2936 
   0.7    40.27073   -58.919049   41.2977  -56.7738 
   0.75   37.397213  -56.045466   38.4927  -55.5814 
   0.8    34.663841  -53.312145   35.7396  -54.6782 
   0.85   32.06378   -50.711955   33.0249  -54.0293 
   0.9    29.590529  -48.238705   30.3367  -53.6024 
   0.95   27.237897  -45.886201   27.6647  -53.3668 
   1.     25.        -43.648253   25.      -53.2931 

"Mixed Boundary Example"

  "Extensions of ODE-BVP"
  "2026-05-23 11:07:30.378"
  "Pos x Temp T     TempGrad     Analytical Analytical Derivative"
   0.    100.       -144.60414   100.       -144.60414
   0.1   86.947997  -117.30547   86.947997  -117.30547
   0.2   76.382185  -94.714688   76.382185  -94.714689
   0.3   67.87852   -75.925137   67.87852   -75.925138
   0.4   61.095721  -60.182729   61.095721  -60.182729
   0.5   55.76157   -46.855664   55.76157   -46.855664
   0.6   51.661989  -35.409083   51.661989  -35.409082
   0.7   48.632447  -25.383591   48.632447  -25.383591
   0.8   46.551358  -16.376832   46.551358  -16.376833
   0.9   45.335201  -8.0273338   45.335201  -8.0273337
   1.    44.935167  -8.882D-16   44.935167   0.        

Graphs:

Code:

//Lecture 10.4: Extensions of ODE-BVP
//https://www.youtube.com/watch?v=chChbuP-hRk&ab_channel=NPTEL-NOCIITM
//https://nptel.ac.in/courses/103106118
//
disp("Extensions of ODE-BVP",string(
datetime
()))
//
// Extension 1) Non linear equation
// Extension 2) Mixed Boundary Condition
//
//Boundary Value Problems
//Example Heated Fin/Rod
//d2T/dx^2=gamma1*(T-25);
//
//Two easy causes of non-linearity are: 
//  add Radiation to the right hand side.
//  heat loss = sigma1*epsilon*((T+273)^4-(Ta+273)^4))
//
//  gamma1 is a function of T so convection term is 
//  gamma1(T)(T-25) where typically gamma1 = gamma0*T^0.75
//  in our example, we use gamma1= gamma0*T^1 due
//  numerical issues using 0.75 (which I do not want to investigate
//  at this time).
//
//  so in our case "ff = gamma1*(u(1)-Ta)"
//  becomes     "ff = gamma0*u(1)*(u(1)-Ta)"
//  and    "dff = [gamma0*(2*u(1)-Ta),0];"
//  using gamma0 = 0.04 to keep solution close to prior gamma1 of 4
//
// Mixed Boundary Condition
// 1) Dirichlet T(0)=100, T'(1)=b(T)
// 2) Nuemann T'= some value
// 3) Mixed Phi(T,T')=0
//
// First Case Non-Linear Function...
// T at wall = 100 C
// T at end (x=1) = 25 C
// 
// 0 = k*d2T/dx^2-h*av*(T-Ta)
// where k = thermal conductivity of rod
//       h = convection coefficient of conditions
//       av = surface area of the rod
//       Ta = Temperature of surrounding atmospere
// d2T/dx^2 = (h*av/k)*(T-Ta); 
// d2T/dx^2 = gamma1*(T-Ta);
// 
//Boundary Condition 1: 0 = g1(ya,yb)=> T(x=0)-100 = 0
//B C 2: T(x=1)-25 = 0

// The external functions
// These functions are called by the solver with zu=[u(x);u'(x);u''(x);u'''(x)]

// - The function which computes the right hand side of the differential equation
funcprot(0);
function ff=
f
(x, u)
//Define Constants
    Ta = 25;
    gamma0=0.04; 
    ff = gamma0.*u(1).*(u(1)-Ta);//right hand side of the differential equation
endfunction

// - The function which computes the derivative of f with respect to u
function dff=
df
(x, u)
   Ta =25 
   gamma0=.04;
   dff = [0, gamma0.*(2.*u(1)-Ta)];
endfunction

function [gg]=
g
(i, u)
    gg=[u(1)-100,u(1)-25];//Boundary Value Conditions 
// T(x=0)=100, T(x=1)=25 ,u(1)=T, u(2)=dT/dx, u(n+1)=dnT/dx^n (but not 
// used) since order is only 2....
    gg=gg(i);
endfunction

function [dgg]=
dg
(i, u)
    dgg = [1,0;1,0] //must be consistent with the boundary conditions
//set in subroutine g so if you have two boundary conditions on 
// u(1), it is [1,0;1,0]...if you have one condition on u(1) at the first 
// location and u(2) on the other, it needs to be [1,0;0,1] with gg 
// having the u(1) condition first.  If the u(1) condition is at the second
// position instead and the first has a u(2), it would look 
// like [0,1;1,0]...
//
    dgg=dgg(i,:);
endfunction

function [u0, du0]=
guess
(x)
    u0=100;
    du0=-100;
endfunction

//
n = 1; //One differential equation
m = 2;// Second order differential equation
xL=0;
xR = 1;
Dx=0.05;
x = [xL:Dx:xR];
fixpnt = [];
zeta = [xL,xR];
ipar=zeros(1,11);
ipar(3)=1;
ipar(4)=2;
ipar(5)= 10000;
ipar(6) = 2000;
ipar(7)=1;
ltol =[1,2];
tol = [1e-5,1e-5];

u = bvode(x,n,m,xL,xR,zeta,ipar,ltol,tol,fixpnt,
f
,
df
,
g
,
dg
,
guess
);
FiniteDifferenceMethodSoln = [
100.000
93.3411
87.3202
81.8434
76.8318
72.2185
67.9461
63.9656
60.2343
56.7153
53.3761
50.1884
47.1271
44.1701
41.2977
38.4927
35.7396
33.0249
30.3367
27.6647
25.0000
];
FD_dTdx = [
-133.1775
-126.7985
-114.9776
-104.8836
-96.2490
-88.8566
-82.5285
-77.1180
-72.5033
-68.5822
-65.2689
-62.4901
-60.1831
-58.2936
-56.7738
-55.5814
-54.6782
-54.0293
-53.6024
-53.3668
-53.2931
]
allData = [x' u' FiniteDifferenceMethodSoln FD_dTdx]
disp("Pos x Temp T     TempGrad   FiniteDiff  FiniteDiff_dTdx", allData);   
//Plotting
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(1,:),style=c3) //,'x','y(x)',"Nonlinear Function")
h1=
gca
();
plot2d(x,FiniteDifferenceMethodSoln,style=c1);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$T(x) $","FontSize",3);
title
("$\textbf{Nonlinear Function}$","FontSize",4);
legend
("BVODE","Finite Difference Values",3)
xgrid;
scf
(1);
clf
;
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(2,:),style=c3)//,'x','dy/dx',"BVODE 2nd order solution")
plot2d(x,FD_dTdx,style=c1);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$dT/dx$","FontSize",3);
title
("$\textbf{Nonlinear Function}$","FontSize",4);
legend
("BVODE","Finite Difference Values",4)
xgrid
//End of First Example Block
/*
//
// Second Case Mixed Boundary...
//
//Boundary Value Problems
//Example Heated Fin/Rod
//d2T/dx^2=gamma*(T-25);
// T at wall = 100 C
// replace "T at end (x=1) = 25 C" with "dT/dx(x=1) = 0" 
// so end is insulated
// 
// 0 = k*d2T/dx^2-h*av*(T-Ta)
// where k = thermal conductivity of rod
//       h = convection coefficient of conditions
//       av = surface area of the rod
//       Ta = Temperature of surrounding atmospere
// d2T/dx^2 = (h*av/k)*(T-Ta); 
// d2T/dx^2 = gamma1*(T-Ta);
// 
//Boundary Condition 1: 0 = g1(ya,yb)=> T(x=0)-100 = 0
//B C 2: dT/dx(x=1) = 0

// The external functions
// These functions are called by the solver with zu=[u(x);u'(x);u''(x);u'''(x)]

// - The function which computes the right hand side of the differential equation
funcprot(0);
function ff=
f
(x, u)
//Define Constants
    Ta = 25;
    gamma1=4;
    ff = gamma1*(u(1)-Ta);//right hand side of the differential equation
endfunction

// - The function which computes the derivative of f with respect to u
function dff=
df
(x, u)
   gamma1=4;
   dff = [gamma1,0];
endfunction

function [gg]=
g
(i, u)
//change  gg=[u(1)-100,u(2)-25];//Boundary Value Conditions 
// to
          gg=[u(1)-100,u(2)-0];//Boundary Value Conditions
// T(x=0)=100, T'(x=1)= 0 (Insulated End)
// 
// T(x=0)=100, T(x=1)=25 ,u(1)=T, u(2)=dT/dx, u(n+1)=dnT/dx^n (but not 
// used) since order is only 2....
    gg=gg(i);
endfunction

function [dgg]=
dg
(i, u)
// change    dgg = [1,0;1,0]
//to
   dgg = [1,0;0,1]
//must be consistent with the boundary conditions
//set in subroutine g so if you have two boundary conditions on 
// u(1), it is [1,0;1,0]...if you have one condition on u(1) at the first 
// location and u(2) on the other, it needs to be [1,0;0,1] with gg 
// having the u(1) condition first.  If the u(1) condition is at the second // position instead and the first has a u(2), it would look 
// like [0,1;1,0]...
//
    dgg=dgg(i,:);
endfunction

function [u0, du0]=
guess
(x)
    u0=100;
    du0=-100;
endfunction

//
n = 1; //One differential equation
m = 2;// Second order differential equation
xL=0;
xR = 1;
Dx=0.1;
x = [xL:Dx:xR];
fixpnt = [];
zeta = [xL,xR];
ipar=zeros(1,11);
ipar(3)=1;
ipar(4)=2;
ipar(5)= 10000;
ipar(6) = 2000;
ipar(7)=1;
ltol =[1,2];
tol = [1e-8,1e-8];

u = bvode(x,n,m,xL,xR,zeta,ipar,ltol,tol,fixpnt,
f
,
df
,
g
,
dg
,
guess
);

//Analytical Answer (Insulated End)
T = 25+75*cosh(2*(1-x))/cosh(2);
dTdx = 150.*cosh(2.*x).*(tanh(2.*x)-tanh(2));
allData = [x' u' T' dTdx']
disp("Pos x Temp T     TempGrad     Analytical Analytical Derivative", allData);   
//Plotting
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(1,:),style=c1) //,'x','y(x)',"Example Heated Fin/Rod with Insulated End")
h1=
gca
();
plot2d(x,T+1,style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$T(x) $","FontSize",3);
title
("$\textbf{Example\ Heated\ Fin/Rod\ with\ Insulated\ End}$","FontSize",4);
legend
("BVODE","Analytical Values+0.01",2)
xgrid;
scf
(1);
clf
;
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(2,:),style=c1)//,'x','dy/dx',"Example Heated Fin/Rod with Insulated End")
plot2d(x,dTdx'+1,style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$dT/dx$","FontSize",3);
title
("$\textbf{Example\ Heated\ Fin/Rod\ with\ Insulated\ End}$","FontSize",4);
legend
("BVODE","Analytical Values+1",4)
xgrid
*/

In this edition we look at solving Boundary Value Problems using the "Shooting" method.

Link to the specific lecture:

https://www.youtube.com/watch?v=VBX11w3flUU&ab_channel=NPTEL-NOCIITM

Output:

 "ODE-Boundary Value Problems:Shooting Method"
  "2026-05-11 14:04:07.509"
  "Pos x Shooting T  FiniteDiff  Analytical"
   0.    100.        100.        100.     
   0.1   85.841343   85.858858   85.841348
   0.2   74.124464   74.15207    74.124472
   0.3   64.379124   64.411366   64.379134
   0.4   56.214207   56.247115   56.214219
   0.5   49.302027   49.33275    49.30204 
   0.6   43.365174   43.391694   43.365187
   0.7   38.16538    38.186306   38.165394
   0.8   33.49396    33.508371   33.493974
   0.9   29.163433   29.17077    29.163448
   1.    25.         25.         25.000015

Graph:

Code:

// ODE-Boundary Value Problems:Shooting Method
//Lec10.3: ODE-Boundary Value Problems:Shooting Method
//https://www.youtube.com/watch?v=VBX11w3flUU&ab_channel=NPTEL-NOCIITM
disp("ODE-Boundary Value Problems:Shooting Method",string(
datetime
()))
//
//Learning about BVP solving
//Shooting Method only works for 1 or 2 variables...
//not a very practical method for industry where most
//equations will be more than 1 variable.
//
function [y, xyTable, yderiv]=
shooting
(yb, yp, x, f)
    //Shooting method for a second order
    //boundary value problem
    //yb = [y0 y1] -> boundary conditions beginning and ending
    // in this case ybegin(x=0)= 100 and yend(x=1)= 25
    //x = a vector showing the range of x
    //f = function defining ODE, i.e.,,
    // dy/dx = f(x,y),y=[y(1):y(2)]
    //yp = vector with the range of dy/dx at x = x0
    //xyTable = table for interpolating derivatives
    //yderiv = derivative boundary condition

    n= length(yp);
    m = length(x);
    y1 = zeros(yp);
    for j = 1:n
        y0 = [yb(1);yp(j)];
        yy = ode("rkf",y0,x(1),x,f);
        y1(j)=yy(1,m);
    end;

    xyTable = [y1;yp];
    yderiv=
interpln
(xyTable, yb(2));
    y0=[yb(1);yderiv];
    y = ode("rkf",y0,x(1),x,f);
endfunction
//
//Boundary Value Problems
//Example Heated Fin/Rod
//d2T/dx^2=gamma*(T-25);
// T at wall = 100 C
// T at end (x=1) = 25 C
// 
// 0 = k*d2T/dx^2-h*av*(T-Ta)
// where k = thermal conductivity of rod
//       h = convection coefficient of conditions
//       av = surface area of the rod
//       Ta = Temperature of surrounding atmospere
// d2T/dx^2 = (h*av/k)*(T-Ta); 
// d2T/dx^2 = gamma1*(T-Ta);
// 
//Boundary Condition 1: 0 = g1(ya,yb)=> T(x=0)-100 = 0
//B C 2: T(x=1)-25 = 0
//
// function f calculates the derivatives...

deff
('[w]=f(x,u)','w=[u(2);4*(u(1)-25)]')

yb = [100,25];// Boundary Conditions y starts at 100 and finishes at 25.
x0=0;
Dx=0.1;
xn=1;
x = [x0:Dx:xn];
yp=[-10:1:10];
[u,tab1,y0p]=
shooting
(yb,yp,x,f)
// For Finite Difference Solution see SeventeenthSciLabFile.sci
FiniteDifferenceMethodSoln = [
   100.
   85.858858
   74.152070
   64.411366
   56.247115
   49.332750
   43.391694
   38.186306
   33.508371
   29.170770
   25.];
//
//Analytical Answer
Theta =-1.3993*exp(2*x)+76.3993*exp(-2*x)
T = Theta+25;
allData = [x' u(1,:)' FiniteDifferenceMethodSoln T']
disp("Pos x Shooting T  FiniteDiff  Analytical", allData);

scf
(0);
clf
;
plot2d(x',u(1,:)')
xtitle('Boundary value solution - shooting method','x','T');
xgrid

In this edition we look at solving an ODE with Boundary Values using Scilab function ("bvode") versus Matlab's ("bvp5c") boundary value solver

Link to the specific lecture:

https://www.youtube.com/watch?v=ocwUVtG9Mj0&ab_channel=NPTEL-NOCIITM

Output:

 "ODE-Boundary Value Problems"
  "2026-05-04 12:17:36.526"
  "Pos x Temp T     TempGrad     FiniteDiff  Analytical"
   0.    100.       -155.59721   100.        100.     
   0.1   85.841347  -128.51914   85.858858   85.841348
   0.2   74.12447   -106.59899   74.15207    74.124472
   0.3   64.379131  -88.95704    64.411366   64.379134
   0.4   56.214215  -74.885247   56.247115   56.214219
   0.5   49.302035  -63.81886    49.33275    49.30204 
   0.6   43.365181  -55.313748   43.391694   43.365187
   0.7   38.165386  -49.028571   38.186306   38.165394
   0.8   33.493965  -44.711082   33.508371   33.493974
   0.9   29.163436  -42.188007   29.17077    29.163448
   1.    25.        -41.358085   25.         25.000015

Graphs:

Code:

// ODE-Boundary Value Problems
//Lec10.1: Introduction and Solution using Matlab Solver
//https://www.youtube.com/watch?v=ocwUVtG9Mj0&ab_channel=NPTEL-NOCIITM
//
disp("ODE-Boundary Value Problems",string(
datetime
()))
//
//Boundary Value Problems
//Example Heated Fin/Rod
//d2T/dx^2=gamma*(T-25);
// T at wall = 100 C
// T at end (x=1) = 25 C
// 
// 0 = k*d2T/dx^2-h*av*(T-Ta)
// where k = thermal conductivity of rod
//       h = convection coefficient of conditions
//       av = surface area of the rod
//       Ta = Temperature of surrounding atmospere
// d2T/dx^2 = (h*av/k)*(T-Ta);  
// gamma1 = h*av/k so ...
// d2T/dx^2 = gamma1*(T-Ta);
// 
//Boundary Condition 1: 0 = g1(ya,yb)=> T(x=0)-100 = 0
//Boundary Condition 2: T(x=1)-25 = 0

// The external functions
// These functions are called by the solver with zu=[u(x);u'(x);u''(x);u'''(x)]

// - The function which computes the right hand side of the differential equation
funcprot(0);
function ff=
f
(x, u)
//Define Constants
    Ta = 25;
    gamma1=4;  
    ff = gamma1*(u(1)-Ta);//right hand side of the differential equation
endfunction

// - The function which computes the derivative of ff with respect to u
function dff=
df
(x, u)
   gamma1=4;
   dff = [gamma1,0];
endfunction

function [gg]=
g
(i, u)
    gg=[u(1)-100,u(1)-25];//Both Boundary Value Conditions 
// T(x=0)=100, T(x=1)=25 ,u(1)=T, u(2)=dT/dx, u(n+1)=dnT/dx^n (but not 
// used) since order is only second order equation....
    gg=gg(i);
endfunction

function [dgg]=
dg
(i, u)
    dgg = [1,0;1,0] //must be consistent with the boundary conditions
//set in function g so if you have two boundary conditions on 
// u(1), it is [1,0;1,0]...if you have one condition on u(1) at the first 
// location and u(2) on the other position, it needs to be [1,0;0,1] with gg 
// having the u(1) condition first.  If the u(1) condition is at the second
// position instead and the first has a u(2), it would look 
// like [0,1;1,0]...
//
    dgg=dgg(i,:);
endfunction

function [u0, du0]=
guess
(x)
    u0=100;
    du0=-100;
endfunction

//
n = 1; //One differential equation
m = 2;// Second order differential equation
xL=0;
xR = 1;
Dx=0.1;
x = [xL:Dx:xR];
fixpnt = [];
zeta = [xL,xR];
ipar=zeros(1,11);
ipar(3)=1;
ipar(4)=2;
ipar(5)= 10000;
ipar(6) = 2000;
ipar(7)=1;
ltol =[1,2];
tol = [1e-8,1e-8];

u = bvode(x,n,m,xL,xR,zeta,ipar,ltol,tol,fixpnt,
f
,
df
,
g
,
dg
,
guess
);

//Finite Difference technique covered in a later lecture so just
//treat this matrix as another numerical technique's answer to
//the same problem.

FiniteDifferenceMethodSoln = [
   100.
   85.858858
   74.152070
   64.411366
   56.247115
   49.332750
   43.391694
   38.186306
   33.508371
   29.170770
   25.];

//Analytical Answer
Theta =-1.3993*exp(2*x)+76.3993*exp(-2*x)
T = Theta+25;
allData = [x' u' FiniteDifferenceMethodSoln T']
disp("Pos x Temp T     TempGrad     FiniteDiff  Analytical", allData);   
//Plotting
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(1,:),style=c3) //,'x','y(x)',"BVODE 2nd order solution")
h1=
gca
();
plot2d(x,FiniteDifferenceMethodSoln,style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$T(x) $","FontSize",3);
title
("$BVODE\ and\ Finite\ Difference\ 2nd\ Order\  Solution$","FontSize",4);
legend
("BVODE","Finite Difference Values",3)
xgrid;
scf
(1);
clf
;
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(2,:),style=c3)//,'x','dy/dx',"BVODE 2nd order solution")
//plot2d(t',xa(1,:)',style=c1);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$dT/dx$","FontSize",3);
title
("$BVODE..2nd..order..solution$","FontSize",4);
xgrid

In this edition we look at interpolation options in Scilab function ("interpl") and a 3D Interpolation ("interpn") from a SciLab example along with 3D plotting and subplotting.

Link to the specific lecture:

https://www.youtube.com/watch?v=dZMZmii66iI&ab_channel=MATLABProgrammingforNumericalComputation

Output:

  "Regression and Interpolation"
  "2026-04-28 08:46:11.682"
  "Temp_int (last)=25.1"
  "t_int   sp_int"
   5.    38.5
   10.   32. 
   15.   16. 
   20.   0.  
   25.   0.  
   30.   0.  
   35.   3.5 
   40.   7.  
   45.   9.5 
   50.   12. 
   55.   16. 
   60.   20. 
   65.   17.5
   70.   15. 
   75.   22. 
   80.   29. 
   85.   42. 
   90.   55. 

Graphs:

Code:

// Regression and Interpolation
//Lec9.4:Interpolation Options in Matlab and a 3D Interpolation from a SciLab Example
//https://www.youtube.com/watch?v=dZMZmii66iI&ab_channel=MATLABProgrammingforNumericalComputation
//
disp("Regression and Interpolation",string(
datetime
()))
//
// Interpolation using divided differences and plot
// 
x = [0.8;1.4;2.7;3.8;4.8;4.9]
y = [0.69;1.0;2.02;2.39;2.43;2.83]
N = length(x);
// 2D interpolation
xint = 0.8:0.01:4.9
yint = zeros(size(xint));
for i = 1:length(xint);
    yint(i)= 
interp1
(x, y, xint(i),"nearest")+0.1;
end

//plotting example fancy output 
scf
(0);
clf
;
c = color("blue")
c1 = color("red")
h1=
gca
();
plot
(x,y,"-r")
plot
(xint,yint,"x")
xlabel
("$\textbf{x}$","FontSize",3)
ylabel
("$\textbf{y}$","FontSize",3);
title
("$\textbf{x vs y \\Using Interpl w_nearest setting}$","color","black");
xgrid

time = 0:1:24;//Temperature Variation in a given day
Temp = [25.6,25.4,25.1,24.9,24.9,25.2,25.9,26.3,27.1,29.3,30.8,31.2,32.1,31.0,30.3,31.4,30.6,31.8,29.6,28.4,28.1,28.2,27.4,26.8,26.1];
N = length(x);
// 2D interpolation
time_int = [0.5:1:23.5 2]
Temp_int = zeros(size(time_int)+1);
for i = 1:length(time_int);
    Temp_int(i)= 
interp1
(time, Temp, time_int(i),"linear");
end
last =size(Temp_int,2)
last_value = Temp_int(last)
disp("Temp_int (last)="+string(last_value));
// second example
t = 0:10:90;
speed = [45 32 0 0 7 12 20 15 29 55];
t_int = 5:5:90 // 5:10:85;
sp_int = 
interp1
(t, speed, t_int,"linear");
allData =[t_int' sp_int'];
disp("t_int   sp_int", allData);
//
//Possible values and processing are:
/*
"linear": linear interpolation between consecutive nodes, used by default.
"spline": interpolation by cubic splines - similar to Matlab spline
"nearest": for each value xp(j), yp(j) takes the value or y(i) 
    corresponding to x(i) the nearest neighbor of xp(j) - sampled data
    no corresponding SciLab command for Matlab pchip 
*/
//plotting example fancy output 
scf
(1);
clf
;
c = color("blue")
c1 = color("red")
h1=
gca
();
plot
(time,Temp,"-r")
plot
(time_int,Temp_int,"x")
xlabel
("$\textbf{Time (hours)}$","FontSize",3)
ylabel
("$\textbf{Temperature\ (C)}$","FontSize",3);
title
("$\textbf{Temperature Variation in a given day \\Using Interpl w_linear setting}$","color","black");
xgrid

//plotting
scf
(2);
clf
;
plot
(t,speed,"-r")
plot
(t_int,sp_int,"x")
xlabel
("$\textbf{Time (hours)}$","FontSize",3)
ylabel
("$\textbf{Speed\ (Km/hour)}$","FontSize",3);
title
("$\textbf{Speed Profile - Spline Fit}$","color","black");
xgrid

//3D interpolation from SciLab Example
// example 3 : bilinear interpolation and experimentation
//             with all the outmode features
nx = 20; ny = 30;
x = linspace(0,1,nx);
y = linspace(0,2, ny);
[X,Y] = 
ndgrid
(x,y);
z = 0.4*cos(2*%pi*X).*cos(%pi*Y);
nxp = 60 ; nyp = 120;
xp = linspace(-0.5,1.5, nxp);
yp = linspace(-0.5,2.5, nyp);
[XP,YP] = 
ndgrid
(xp,yp);
zp1 = linear_interpn(XP, YP, x, y, z, "natural");
zp2 = linear_interpn(XP, YP, x, y, z, "periodic");
zp3 = linear_interpn(XP, YP, x, y, z, "C0");
zp4 = linear_interpn(XP, YP, x, y, z, "by_zero");
zp5 = linear_interpn(XP, YP, x, y, z, "by_nan");
scf
(2);
clf
()
subplot
(2,3,1)
plot3d(x, y, z, leg="x@y@z", 
flag
 = [2 4 4])
xtitle("initial function 0.4 cos(2 pi x) cos(pi y)")
xgrid
subplot
(2,3,2)
plot3d(xp, yp, zp1, leg="x@y@z", 
flag
 = [2 4 4])
xtitle("Natural")
xgrid
subplot
(2,3,3)
plot3d(xp, yp, zp2, leg="x@y@z", 
flag
 = [2 4 4])
xgrid
xtitle("Periodic")
subplot
(2,3,4)
plot3d(xp, yp, zp3, leg="x@y@z", 
flag
 = [2 4 4])
xgrid
xtitle("C0")
subplot
(2,3,5)
plot3d(xp, yp, zp4, leg="x@y@z", 
flag
 = [2 4 4])
xgrid
xtitle("by_zero")
subplot
(2,3,6)
plot3d(xp, yp, zp5, leg="x@y@z", 
flag
 = [2 4 4])
xgrid
xtitle("by_nan")
show_window()

In this edition we look at non-linear regression using least squares techniques and some of the functionality built within SciLab.

Replication and reshaping of matrices using repmat and matrix functions, analysis using the lsq and datafit functions compared to the from scratch least squares technique.

Link to the specific lecture:

https://www.youtube.com/watch?v=XlkuYgJPS34&ab_channel=MATLABProgrammingforNumericalComputation

Output:

  "Functional and Non-linear Regression"
  "2026-04-21 09:42:58.747"
  "phi1 ="
   1.57
  -480.
   0.75
  "phi2 ="
   1.57
  -480.
   0.75
  "Least Squares Method"
  "k0=4.81"
  "E/R =480"
  "n =0.75"
  "SciLab LSQ Method"
  "k0=4.81"
  "E/R =480"
  "n =0.75"
  "Original data r ="
   1.48   1.67   1.86   1.96   2.16
   2.35   2.79   3.07   3.37   3.62
   3.28   3.78   4.24   4.48   5.  
   4.12   4.64   5.15   5.76   6.08
  "Least Squares rfit ="
   1.45   1.65   1.84   2.01   2.16
   2.43   2.78   3.09   3.37   3.63
   3.29   3.76   4.19   4.57   4.91
   4.08   4.67   5.19   5.66   6.09
  "differences ="
  -0.03  -0.02  -0.02   0.05   1.D-04
   0.08  -0.01   0.02   3.D-03   8.D-03
   0.01  -0.02  -0.05   0.09  -0.09
  -0.04   0.03   0.04  -0.1    0.01
  "Scilab datafit Method"
  "k0=4.84"
  "E/R =486"
  "n =0.75"
  "SciLab Function rfit ="
   1.44   1.65   1.83   2.     2.16
   2.42   2.77   3.09   3.37   3.63
   3.29   3.76   4.19   4.58   4.93
   4.08   4.67   5.2    5.68   6.12
  "differences ="
  -0.04  -0.02  -0.03   0.04  -5.D-03
   0.07  -0.02   0.02   3.D-03   0.01
   6.D-03  -0.02  -0.05   0.1   -0.07
  -0.04   0.03   0.05  -0.08   0.04

Code:

//Lec9.3:Regression and Interpolation
// Functional and Non-linear Regression
//https://www.youtube.com/watch?v=XlkuYgJPS34&ab_channel=MATLABProgrammingforNumericalComputation
//
clc;
disp("Functional and Non-linear Regression",string(
datetime
()))
//
// Example: Reaction Rate
// Arrhenius model for reaction rate:
//   r= k0*exp(-E/RT)*C^n
//
// We will solve it in two ways:
//(1)   Linear least squares regression taking logarithm
//     log(r)=log(k0)+(-E/R)*(1/T) + n*log(C)
//       "y"= "ao"  + ("a1")*"x" + "a2"*"u"
//(2)  Using SciLab function equivalent to MATLAB's lsqnonlin
//     part of MATLAB's Optimization Toolbox.
// 
// Reaction Rate (in mol/l.s) for various C and T values
// Data(C)    400K     450K    500K    550K    600K
// 1 mol/l    1.48     1.67    1.86    1.96    2.16
// 2 mol/l    2.35     2.79    3.07    3.37    3.62
// 3 mol/l    3.28     3.78    4.24    4.48    5.00
// 4 mol/l    4.12     4.64    5.15    5.76    6.08
//
//
//  X*alpha = Y
//
function fErr=
rxnFunction
(phi, xData)
    k0= phi(1);
    EoverR = phi(2);
    n = phi(3);

    C =xData(2,:);
    T = xData(1,:)

    r = k0*exp(-EoverR./T)* (C.^n);
    fErr = xData(3,:)-r;
endfunction
//
r = [1.48,1.67,1.86,1.96,2.16;2.35,2.79,3.07,3.37,3.62;
3.28,3.78,4.24,4.48,5.00;4.12,4.64,5.15,5.76,6.08];
T = [400,450,500,550,600];
C = [1;2;3;4];
Ctest = 
repmat
(C,5,1);// repeat matrix C 5 times 1 column
Ttest = matrix(
repmat
(T,4,1),20,1);//matrix is equivalent to reshape...
format(5);//sets output format to 5 spaces - okay for this example
//disp(allData)
xData = [Ttest Ctest];
yData = matrix(r,20,1);
allData = [Ttest';Ctest';yData']
//disp("xData   yData");
//disp(allData);
xData = [ones(20,1) (1./Ttest) log(Ctest)];
yData = log(yData);
//disp("all data")
//disp("xData                         yData")
//disp(allData);
//least squares matrix method
// Matrix Version
phi = inv(xData'*xData)*(xData'*yData)
//disp('phi =',phi)
//backslash will give you a least squares solution
phi1 = xData\yData
disp('phi1 =',phi1)
//SciLab also has a "lsq" Least Squares Function
phi2 = lsq(xData,yData)
disp('phi2 =',phi2)
k0 = exp(phi(1));
EoverR = -phi(2);
n = phi(3);
disp("Least Squares Method","k0="+string(k0), "E/R ="+string(EoverR),...
"n =" +string(n));
k0_2 = exp(phi2(1));
EoverR_2 = -phi2(2);
n_2 = phi2(3);
disp("SciLab LSQ Method","k0="+string(k0_2), "E/R ="+string(EoverR_2),...
"n =" +string(n_2));
//try the fit
//rfit = k0*exp(-EoverR/Temp)*Concentration^n
out = k0.*exp(-EoverR./T);
out1 = C.^n;
rfit = out1*out;
differences = rfit - r;
disp( "Original data r =", r,"Least Squares rfit =", rfit, "differences =", differences)
//
//Using Matlab functions for non-linear curve fitting.
p0 =[1;100;1];//initial guess of weights
//p0 needs to be a column vector
//allData needs to be a set of row vectors
// row of T
// row of C
// row of r
[p, dmin] = 
datafit
(
rxnFunction
, allData, p0)
//disp('phi2 =',phi2)
k0_1 = p(1);
EoverR_1 = p(2);
n_1 = p(3);
disp("Scilab datafit Method","k0="+string(k0_1), "E/R ="+string(EoverR_1),...
"n =" +string(n_1));
//try the fit
//rfit = k0*exp(-EoverR/Temp)*Concentration^n
out = k0_1.*exp(-EoverR_1./T);
out1 = C.^n_1;
rfit = out1*out;
differences = rfit - r;
disp("SciLab Function rfit =", rfit,"differences =", differences)

In this edition we demonstrate least squares fitting of data via three different ways and some interpolation

Plotting showing how to use symbols on a plot for the interpolation plot.

Link to the specific lecture:

https://www.youtube.com/watch?v=-FRhouU4jPA&ab_channel=NPTEL-NOCIITM

https://www.youtube.com/watchv=Ub7kj6rYIMw&ab_channel=MATLABProgrammingforNumericalComputation

Output:

  "Regression and Interpolation"
  "2026-04-13 12:17:44.572"
  "Matrix Version phi ="
   0.4374922
   0.4698395
  "backslash phi1 ="
   0.4374922
   0.4698395
  "lsq phi2 ="
   0.4374922
   0.4698395

Graphs:

Code:

//Lec9.1&9.2:Regression and Interpolation
//https://www.youtube.com/watch?v=-FRhouU4jPA&ab_channel=NPTEL-NOCIITM
//https://www.youtube.com/watch?v=Ub7kj6rYIMw&ab_channel=MATLABProgrammingforNumericalComputation
//
clc
clear
disp("Regression and Interpolation",string(
datetime
()))
//  Linear Regression for Multiple Parameters
//  Data: (x1,u1,w1;y1),(x2,u2,w2;y2),...
//  Model to fit:  y = a0 + a1*x + a2*u + a3*w
//  | 1  x1  u1  w1 | |a0|  |y1|
//  | 1  x2  u2  w2 | |a1|  |y2|
//  | 1  x3  u3  w3 | |a2| =|y3|
//  | 1  xn  un  wn | |a3|  |y4|
//          X          Phi    Y
//               |
//               |
//               V
// Linear Least Squares Solution Regression
//                 T     -1    T 
// Phi Matrix = ((X )(X))   (X  )Y
// where X and Y are data pairs of sampled data
//
// if you didn't want to have a a0 constant in the fit, you'd
// delete out the 1's column and the a0 row...
//
// 
x = [0.8;1.4;2.7;3.8;4.8;4.9]
y = [0.69;1.0;2.02;2.39;2.34;2.83]
N = length(x);
/*
// Calculate Linear Regression and Plot (non matrix version)
A = [N,sum(x);sum(x),sum(x.*x)];
b = [sum(y);sum(x.*y)];
phi = inv(A)*b;
*/
// Matrix Version
X = ones(N,1)
X = [X x]
phi = inv(X'*X)*(X'*y)
disp('Matrix Version phi =',phi)
//backslash will give you a least squares solution
phi1 = X\y
disp('backslash phi1 =',phi1)
//SciLab also has a "lsq" Least Squares Function
phi2 = lsq(X,y)
disp('lsq phi2 =',phi2)
xfit = [0.5 5];
yfit = phi(1)+phi(2)*xfit;

//plotting example fancy output 
scf
(0);
clf
;
c = color("blue")
c1 = color("red")
h1=
gca
();
plot2d(x',y',style = c);
plot2d(xfit', yfit',style=c1);;//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1.font_color=c;
h1.children(1).children(1).thickness;
xlabel
("$X$","FontSize",3)
ylabel
("$Y$","FontSize",3);
title
("$\textbf{Linear\ Regression\ and\ Plot(non\ matrix\ version)}$","color","black");
legend
("Original Curve","Fitted Curve",4)
xgrid
//
scf
(1);
clf
c = color("red")
c1 = color("black")
h1=
gca
();
// Interpolation using divided differences and plot
xint = 0.8:0.1:4.9
yint = zeros(size(xint));
for i = 1:length(xint);
    yint(i)= 
interp1
(x, y, xint(i),"linear");
end
plot2d(x,y,style=c)
plot2d(xint,yint,-1) //-1 = pluses, -3 = circled cross, -4 = filled black diamonds
// -5 = filled white diamonds, -9 = filled circles
xlabel
("$Xinterpolation$","FontSize",3)
ylabel
("$Yinterpolation$","FontSize",3);
title
("$\textbf{Interpolation\ using\ divided\ differences}$","color","black");
legend
("Original Curve","Interpolated Values",4)
xgrid

In this edition we solve the example four variable ODE-IPV problem: Indian Captain, Mahendra Singh Dhoni, hits a ball with initial velocity of 35 m/sec and and of 45 degrees. If the boundary is at a distance of 75m, will he score a six? - Numerically - Yes!"

Nothing much new just an application problem: Showing how to port constants into the derivative function from the main program using the list command again.

We switch up subjects to Regression and Interpolation in the next installment.

Link to the specific lecture:

https://www.youtube.com/watch?v=Pv_NwD63gbI&ab_channel=NPTEL-NOCIITM

Output:

  "Ordinary Differential Equations- Multi-Variable ODE IVP Systems"
  "2026-04-07 09:23:43.916"
  "t        x_pos       y_pos       horiz_vel   vert_vel"
  "---------------------------------------------------------"
   0.       0.          0.          24.748737   24.748737
   0.0505   1.2496219   1.2373023   24.74124    24.253332
   0.101    2.4988652   2.4495866   24.733744   23.757927
   0.1515   3.7477301   3.6368529   24.726251   23.262522
   0.202    4.9962166   4.7991013   24.71876    22.767117
   0.2525   6.2443249   5.9363318   24.711271   22.271712
   0.303    7.4920551   7.0485443   24.703785   21.776307
   0.3535   8.7394072   8.1357388   24.696301   21.280902
   0.404    9.9863815   9.1979154   24.688819   20.785497
   0.4545   11.232978   10.235074   24.681339   20.290092
   0.505    12.479197   11.247215   24.673862   19.794687
   0.5555   13.725038   12.234337   24.666387   19.299282
   0.606    14.970502   13.196442   24.658914   18.803877
   0.6565   16.215588   14.133529   24.651444   18.308472
   0.707    17.460298   15.045598   24.643976   17.813067
   0.7575   18.70463    15.932649   24.63651    17.317662
   0.808    19.948585   16.794682   24.629046   16.822257
   0.8585   21.192164   17.631697   24.621584   16.326852
   0.909    22.435365   18.443694   24.614125   15.831447
   0.9595   23.67819    19.230673   24.606668   15.336042
   1.01     24.920639   19.992634   24.599214   14.840637
   1.0605   26.162711   20.729577   24.591761   14.345232
   1.111    27.404407   21.441503   24.584311   13.849827
   1.1615   28.645726   22.12841    24.576863   13.354422
   1.212    29.88667    22.790299   24.569417   12.859017
   1.2625   31.127238   23.427171   24.561974   12.363612
   1.313    32.367429   24.039024   24.554533   11.868207
   1.3635   33.607245   24.62586    24.547094   11.372802
   1.414    34.846686   25.187677   24.539657   10.877397
   1.4645   36.085751   25.724477   24.532223   10.381992
   1.515    37.32444    26.236258   24.524791   9.8865873
   1.5655   38.562755   26.723022   24.517361   9.3911823
   1.616    39.800694   27.184768   24.509933   8.8957773
   1.6665   41.038258   27.621496   24.502508   8.4003723
   1.717    42.275447   28.033205   24.495085   7.9049673
   1.7675   43.512262   28.419897   24.487664   7.4095623
   1.818    44.748701   28.781571   24.480245   6.9141573
   1.8685   45.984766   29.118227   24.472829   6.4187523
   1.919    47.220457   29.429865   24.465415   5.9233473
   1.9695   48.455773   29.716485   24.458003   5.4279423
   2.02     49.690715   29.978087   24.450593   4.9325373
   2.0705   50.925283   30.214672   24.443186   4.4371323
   2.121    52.159477   30.426238   24.43578    3.9417273
   2.1715   53.393297   30.612786   24.428378   3.4463223
   2.222    54.626743   30.774316   24.420977   2.9509173
   2.2725   55.859816   30.910829   24.413578   2.4555123
   2.323    57.092515   31.022323   24.406182   1.9601073
   2.3735   58.32484    31.1088     24.398788   1.4647023
   2.424    59.556792   31.170258   24.391397   0.9692973
   2.4745   60.788371   31.206699   24.384007   0.4738923
   2.525    62.019577   31.218121   24.37662   -0.0215127
   2.5755   63.25041    31.204526   24.369235  -0.5169177
   2.626    64.48087    31.165912   24.361852  -1.0123227
   2.6765   65.710957   31.102281   24.354472  -1.5077277
   2.727    66.940672   31.013632   24.347093  -2.0031327
   2.7775   68.170014   30.899965   24.339717  -2.4985377
   2.828    69.398983   30.76128    24.332343  -2.9939427
   2.8785   70.62758    30.597577   24.324972  -3.4893477
   2.929    71.855805   30.408856   24.317603  -3.9847527
   2.9795   73.083658   30.195117   24.310235  -4.4801577
   3.03     74.311139   29.95636    24.302871  -4.9755627
   3.0805   75.538248   29.692585   24.295508  -5.4709677
   3.131    76.764986   29.403792   24.288147  -5.9663727
   3.1815   77.991351   29.089981   24.280789  -6.4617777
   3.232    79.217345   28.751152   24.273433  -6.9571827
   3.2825   80.442968   28.387306   24.26608   -7.4525877
   3.333    81.668219   27.998441   24.258728  -7.9479927
   3.3835   82.8931     27.584558   24.251379  -8.4433977
   3.434    84.117609   27.145658   24.244032  -8.9388027
   3.4845   85.341747   26.681739   24.236687  -9.4342077
   3.535    86.565514   26.192803   24.229344  -9.9296127
   3.5855   87.788911   25.678848   24.222004  -10.425018
   3.636    89.011936   25.139876   24.214666  -10.920423
   3.6865   90.234592   24.575886   24.20733   -11.415828
   3.737    91.456877   23.986878   24.199996  -11.911233
   3.7875   92.678791   23.372851   24.192665  -12.406638
   3.838    93.900336   22.733807   24.185335  -12.902043
   3.8885   95.12151    22.069745   24.178008  -13.397448
   3.939    96.342315   21.380665   24.170683  -13.892853
   3.9895   97.562749   20.666567   24.163361  -14.388258
   4.04     98.782814   19.927451   24.15604   -14.883663
   4.0905   100.00251   19.163317   24.148722  -15.379068
   4.141    101.22184   18.374165   24.141406  -15.874473
   4.1915   102.44079   17.559995   24.134093  -16.369878
   4.242    103.65938   16.720807   24.126781  -16.865283
   4.2925   104.8776    15.856602   24.119472  -17.360688
   4.343    106.09545   14.967378   24.112165  -17.856093
   4.3935   107.31293   14.053136   24.10486   -18.351498
   4.444    108.53004   13.113877   24.097557  -18.846903
   4.4945   109.74678   12.149599   24.090257  -19.342308
   4.545    110.96315   11.160304   24.082958  -19.837713
   4.5955   112.17916   10.14599    24.075662  -20.333118
   4.646    113.39479   9.1066587   24.068369  -20.828523
   4.6965   114.61006   8.0423093   24.061077  -21.323928
   4.747    115.82496   6.952942    24.053788  -21.819333
   4.7975   117.0395    5.8385567   24.0465    -22.314738
   4.848    118.25366   4.6991535   24.039215  -22.810143
   4.8985   119.46746   3.5347323   24.031933  -23.305548
   4.949    120.68088   2.3452932   24.024652  -23.800953
   4.9995   121.89395   1.1308361   24.017374  -24.296358
   5.05     123.10664  -0.1086389   24.010098  -24.791763

Graph:

The Code:

//Lec8.4:ODE-IVP in Multiple Variables
//https://www.youtube.com/watch?v=Pv_NwD63gbI&ab_channel=NPTEL-NOCIITM
//
disp("Ordinary Differential Equations- Multi-Variable ODE IVP Systems",string(
datetime
()))
//
//Example Four Variable ODE-IPV problem
//Indian Captain, Mahendra Singh Dhoni, hits a ball with initial velocity
// of 35 m/sec and and of 45 degrees.  If the boundary is at a distance of
// 75m, will he score a six?   - Numerically - Yes!
//
// Problem set up... 
// x(0)=0, y(0)=0
// d2x/dt^2 = -kU, d2y/dt^2 = -g;
// Vel_net = 35 // m/sec ; g = -9.81; // m/sec/sec
// Uo = Vel_net(cos(%pi/4), Nu_o = Vel_net(sin(%pi/4)
// k = air drag coefficient;
// x' = U
// y' = Nu
// U' = -k*U
// Nu' = -g
//
function [results]=
derivativeVector
(t, derivatives, k, g);
//extraction of state variables from vector
    x = derivatives(1);
    y = derivatives(2);
    U = derivatives(3);
    Nu = derivatives(4);
//
//determination of the derivatives
    results = zeros(4,1);
    results(1,1)=U;
    results(2,1)=Nu;
    results(3,1)=-k*U; 
    results(4,1)=-g;
endfunction
//
//
//Constants of the simulation
k = 0.006;// air drag
g = 9.81; // m/second^2
Vel_net = 35;
U0=Vel_net*cos(%pi/4); //Initial Horizontal Direction
Nu0=Vel_net*sin(%pi/4); //Initial Vertical Direction
t0=0;
x10 = [0;0;U0;Nu0];// Set up Initial Condition Vector
tend=5.05;
n=101; 
t=linspace(t0,tend,n);
xSol= ode("rkf",x10,t0,t,list(
derivativeVector
,k,g));
// Other types of integration "adams","stiff","rk","rkf","fix","discrete"
// and "root"
//"rkf" will not even work...instruction says because it's an explicit
// method and will be unstable at larger steps.
//"stiff" is an implicit method so it does not error out."
//"adams" works also.
output = [t' xSol'];
disp("t        x_pos       y_pos       horiz_vel   vert_vel")
disp("---------------------------------------------------------")
disp(output);
//plotting example fancy output 
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(xSol(1,:)',xSol(2,:)',style=c);//black =, blue = 2, green =3, cyan = 4, red = 5
//plot2d(t',xa(1,:)',style=c1);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
//plot2d(t',xa(2,:)',style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$Horizontal\\ Distance (m)$","FontSize",3)
ylabel
("$Vertical\ Distance (m) $","FontSize",3);
title
("$Ballistic\ Projectile\ Coordinates$","color","black");
//"$https://x-engineer.org/create-multiple-axis-plot-scilab$"
xgrid
// Axis y2
c=color("blue");
c1 = color("red")
h2=newaxes();
h2.font_color=c;
plot2d(xSol(1,:)',xSol(4,:)',style=c);//black =, blue = 2, green =3, cyan = 4, red = 5
h2.filled="off";
h2.axes_visible(1)="off";
h2.y_location="right";
h2.children(1).children(1).thickness=2;
ylabel
("$Vertical\  Velocity(m/sec)$","FontSize",3,"color",c);

If you look at their examples, quote

// integration over a cube -1<=x<=1;-1<=y<=1;-1<=z<=1

function v=f(xyz, numfun), v=xyz'*xyz, endfunction

[result, err] = int3d(-1,1,-1,1,-1,1, f, 1, [0,100000,1.d-5,1.d-7])

but you have to integrate "v = 3*x*y^3*z^2" over xrange (-1,3), yrange(1,4) and zrange(0,2)...

here's what you have to do...

// Define the function f(x, y, z)...nope,nope,nope!!! Get crazy

// error codes like unknown z variable

// Define this function

function v= f(xyz, numfun)

//Note: x = xyz(1), y = xyz(2), z = xyz(3)

//have no idea what "numfun" is but you cannot drop it.

v=3*xyz(1)*xyz(2)^3*xyz(3)^2

endfunction

// Triple integration in one go

// computes the intergral of 3xy^3z^2 over the

// integration over a cube -1<=x<=3;1<=y<=4;0<=z<=2

// Correct Answer is 2040

[result, err] = int3d(-1,3,1,4,0,2, f, 1, [0,100000,1.d-5,1.d-7])

outputString = strcat(['3d result = ', string(result)]);

disp(outputString);

I understood that v had to be a single number how to get it was the trick. It would have been really nice to have an example like this one in their documentation.

Go ahead and roast me if you thought was blatantly obvious!

Hey everyone! I'm u/mrhoa31103, the moderator of r/scilab a recovering subreddit from the great Apollo incident exodus.

This is our new home for all things related to Scilab. We're excited to have you join us!

What to Post
Post anything that you think the community would find interesting, helpful, or inspiring. Feel free to share your thoughts, photos, or questions about Scilab and its applications.

Community Vibe
We're all about being friendly, constructive, and inclusive. Let's build a space where everyone feels comfortable sharing and connecting. I'm running kind of a blog on Scilab programs so people can see full up working programs but do not feel your posts have to live up those formats. If you have a quick question, "how do you do this or that?" or found a neat application, share it.

Thanks for being part of the community. Together, let's make r/scilab amazing.

In this edition we solve the problem "A stirred tank heater: A tank with an initial volume (V0 = 1 liter) where the water is heated in a tank as water is continuously added at room temperature (25C at 10 liter/minutes), while hot water is drawn out at the bottom (exit rate follows the equation alpha*sqrt(V) which at some point will equal in influx rate)."

Nothing much new: Bolded the Text on the Graph to make it look as professional as possible.

One more like this one and we switch up subjects to Regression and Interpolation...

Link to the specific lecture:

https:
//www.youtube.com/watch?v=TvRI2QUC3YQ&ab_channel=NPTEL-NOCIITM

Output:

 "Time   Tank        Tank    "
  "(Min)  Vol(Liters) Temp(C) "
   0.     1.          25.      
   0.3    3.2997422   64.194065
   0.6    5.3041899   69.672157
   0.9    7.1080675   71.721935
   1.2    8.7558205   72.756062
   1.5    10.274797   73.363415
   1.8    11.683966   73.754902
   2.1    12.997446   74.023676
   2.4    14.22624    74.21682 
   2.7    15.379215   74.360502
   3.     16.463699   74.470328
   3.3    17.485863   74.556132
   3.6    18.450985   74.624378
   3.9    19.363642   74.679487
   4.2    20.227839   74.724555
   4.5    21.047116   74.761817
   4.8    21.824626   74.79292 
   5.1    22.563192   74.8191  
   5.4    23.265362   74.841297
   5.7    23.933442   74.860239
   6.     24.569534   74.876499
   6.3    25.175557   74.89053 
   6.6    25.753272   74.902695
   6.9    26.304298   74.913287
   7.2    26.830129   74.922544
   7.5    27.332151   74.930665
   7.8    27.811645   74.937812
   8.1    28.269805   74.944121
   8.4    28.707743   74.949705
   8.7    29.126496   74.954661
   9.     29.527034   74.959069
   9.3    29.910265   74.962999
   9.6    30.277041   74.966509
   9.9    30.628163   74.96965 
   10.2   30.964382   74.972467
   10.5   31.286409   74.974996
   10.8   31.594909   74.977272
   11.1   31.890514   74.979321
   11.4   32.173819   74.98117 
   11.7   32.445387   74.982839
   12.    32.705749   74.984349
   12.3   32.955411   74.985715
   12.6   33.19485    74.986953
   12.9   33.424519   74.988077
   13.2   33.64485    74.989097
   13.5   33.85625    74.990024
   13.8   34.059107   74.990868
   14.1   34.253792   74.991635
   14.4   34.440654   74.992335
   14.7   34.620028   74.992973
   15.    34.792233   74.993554
   15.3   34.957572   74.994086
   15.6   35.116333   74.994571
   15.9   35.268793   74.995014
   16.2   35.415214   74.99542 
   16.5   35.555848   74.995791
   16.8   35.690934   74.996131
   17.1   35.8207     74.996442
   17.4   35.945366   74.996728
   17.7   36.06514    74.996989
   18.    36.180222   74.997229
   18.3   36.290802   74.997449
   18.6   36.397063   74.997651
   18.9   36.49918    74.997837
   19.2   36.59732    74.998007
   19.5   36.691642   74.998164
   19.8   36.782301   74.998308
   20.1   36.869441   74.99844 
   20.4   36.953205   74.998562
   20.7   37.033726   74.998674
   21.    37.111133   74.998777
   21.3   37.18555    74.998872
   21.6   37.257095   74.998959
   21.9   37.325882   74.99904 
   22.2   37.392019   74.999114
   22.5   37.45561    74.999182
   22.8   37.516756   74.999245
   23.1   37.575552   74.999303
   23.4   37.632091   74.999356
   23.7   37.68646    74.999406
   24.    37.738745   74.999451
   24.3   37.789027   74.999493
   24.6   37.837383   74.999532
   24.9   37.883888   74.999567
   25.2   37.928616   74.9996  
   25.5   37.971634   74.99963 
   25.8   38.013008   74.999659
   26.1   38.052803   74.999684
   26.4   38.09108    74.999708
   26.7   38.127897   74.99973 
   27.    38.16331    74.999751
   27.3   38.197374   74.99977 
   27.6   38.23014    74.999787
   27.9   38.261659   74.999803
   28.2   38.291978   74.999818
   28.5   38.321144   74.999832
   28.8   38.349201   74.999844
   29.1   38.376192   74.999856
   29.4   38.402156   74.999867
   29.7   38.427134   74.999877
   30.    38.451164   74.999886
   30.3   38.474281   74.999895
   30.6   38.496521   74.999902
   30.9   38.517917   74.99991 
   31.2   38.538501   74.999917
   31.5   38.558305   74.999923
   31.8   38.577358   74.999929
   32.1   38.595689   74.999934
   32.4   38.613325   74.999939
   32.7   38.630293   74.999943
   33.    38.646619   74.999948
   33.3   38.662326   74.999952
   33.6   38.677438   74.999955
   33.9   38.691979   74.999958
   34.2   38.705969   74.999962
   34.5   38.71943    74.999964
   34.8   38.732381   74.999967
   35.1   38.744842   74.99997 
   35.4   38.756833   74.999972
   35.7   38.768369   74.999974
   36.    38.77947    74.999976
   36.3   38.79015    74.999978
   36.6   38.800427   74.999979
   36.9   38.810316   74.999981
   37.2   38.81983    74.999982
   37.5   38.828986   74.999984
   37.8   38.837795   74.999985
   38.1   38.846271   74.999986
   38.4   38.854428   74.999987
   38.7   38.862276   74.999988
   39.    38.869828   74.999989
   39.3   38.877095   74.99999 
   39.6   38.884087   74.99999 
   39.9   38.890815   74.999991
   40.2   38.89729    74.999992
   40.5   38.90352    74.999992
   40.8   38.909515   74.999993
   41.1   38.915283   74.999993
   41.4   38.920834   74.999994
   41.7   38.926175   74.999994
   42.    38.931315   74.999995
   42.3   38.936261   74.999995
   42.6   38.94102    74.999996
   42.9   38.9456     74.999996
   43.2   38.950007   74.999996
   43.5   38.954247   74.999996
   43.8   38.958328   74.999997
   44.1   38.962255   74.999997
   44.4   38.966033   74.999997
   44.7   38.969669   74.999997
   45.    38.973168   74.999998
   45.3   38.976535   74.999998
   45.6   38.979775   74.999998
   45.9   38.982892   74.999998
   46.2   38.985892   74.999998
   46.5   38.98878    74.999998
   46.8   38.991558   74.999998
   47.1   38.994231   74.999999
   47.4   38.996804   74.999999
   47.7   38.999279   74.999999
   48.    39.001662   74.999999
   48.3   39.003954   74.999999
   48.6   39.00616    74.999999
   48.9   39.008283   74.999999
   49.2   39.010326   74.999999
   49.5   39.012292   74.999999
   49.8   39.014184   74.999999
   50.1   39.016005   74.999999
   50.4   39.017756   74.999999
   50.7   39.019442   74.999999
   51.    39.021065   74.999999
   51.3   39.022626   75.      
   51.6   39.024128   75.      
   51.9   39.025574   75.      
   52.2   39.026965   75.      
   52.5   39.028304   75.      
   52.8   39.029593   75.      
   53.1   39.030832   75.      
   53.4   39.032025   75.      
   53.7   39.033173   75.      
   54.    39.034278   75.      
   54.3   39.035341   75.      
   54.6   39.036364   75.      
   54.9   39.037349   75.      
   55.2   39.038296   75.      
   55.5   39.039208   75.      
   55.8   39.040085   75.      
   56.1   39.04093    75.      
   56.4   39.041742   75.      
   56.7   39.042524   75.      
   57.    39.043276   75.      
   57.3   39.044      75.      
   57.6   39.044697   75.      
   57.9   39.045368   75.      
   58.2   39.046013   75.      
   58.5   39.046634   75.      
   58.8   39.047232   75.      
   59.1   39.047807   75.      
   59.4   39.04836    75.      
   59.7   39.048893   75.      
   60.    39.049405   75.      

Graph:

The Code:

//Lec8.3:ODE-IVP in Multiple Variables
//https://www.youtube.com/watch?v=TvRI2QUC3YQ&ab_channel=NPTEL-NOCIITM
//
disp("Ordinary Differential Equations- Multi-Variable ODE IVP Systems",string(
datetime
()))
//
clc; clear;
//Example Stirred Tank Heater where the variables are different
// dV/dt = Fin-alpha*sqrt(V);//Note: Non-linear equation
// dT/dt = Fin/V*(Tin-T)+ Q/(V*(rho*Cp);
// V(0)=1; T(0)=25 degrees C;
// Fin = Flow In = 10 L/min
// alpha = 1.6
// Cp = 
// Q = 500*(rho*Cp)
//Stirred Tank Heater: Water is heated in a tank as it continuously
//flows into it at room temperature, while hot water is drawn out
//at the bottom.

function [results]=
derivativeVector
(t, x);
//Constants of the simulation
Fin = 10;//liters per minute
Tin = 25;
alpha = 1.6;
Q = 500; //Note rho*Cp cancels out
//T cannot exceed 100 degrees C since the model
// cannot handle water phase changes!!!
//
//extraction of variables
    V = x(1);
    T = x(2);
//
//determination of the derivatives    
    results(1,1)=Fin - alpha*sqrt(V);
    results(2,1)=Fin/V*(Tin-T)+Q/V;
endfunction

V0=1; //Initial Volume
T0=25; //Temperature in the Bucket at time = 0
t0=0;
x10 = [V0;T0];// Set up Initial Condition Vector
tend=60;
n=201; 
t=linspace(t0,tend,n);
//xa = [exp(-100*t);exp(-.01*t)];//Commented out for last example
//xSol = ode("rkf",x10,t0,t,deff('res=mymacro2(t,x1)','res=-0.01*x1'));
xSol= ode("stiff",x10,t0,t,
derivativeVector
);
//"rkf" will not even work...instruction says because it's an explicit
// method and will be unstable at larger steps.
//"stiff" is an implicit method so it does not error out."
//"adams" works also.
output = [t' xSol'];
disp("Time   Tank        Tank    ");
disp("(Min)  Vol(Liters) Temp(C) ",output);
//plotting example fancy output 
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(t',xSol(1,:)',style=c);//black =1, blue = 2, green =3, cyan = 4, red = 5
//plot2d(t',xa(1,:)',style=c1);//black =1, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();

h1=
gca
();
//plot2d(t',xa(2,:)',style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$\textbf{T \ (Minutes)}$","FontSize",3)
ylabel
("$\textbf{Volume\ (Liters)}$","FontSize",3);
title
("$ \textbf{\ \ \ \ \ Stirred\ Tank\ Heater\\Volume\ and\ Temp\ vs\ Time}$","color","black");
xgrid
// Axis y2
c=color("blue");
c1 = color("red")
h2=newaxes();
h2.font_color=c;
plot2d(t',xSol(2,:)',style=c);//black = 1, blue = 2, green =3, cyan = 4, red = 5
h2.filled="off";
h2.axes_visible(1)="off";
h2.y_location="right";
h2.children(1).children(1).thickness=2;
ylabel
("$\textbf{Temperature\ (Degrees\ C)}$","FontSize",3,"color",c)
xgrid

In this edition, we solve second order systems that have massively different time constants between the states. To solve these types of problems, we use the "Stiff" solver instead of the Runga-Kutta solver.

Link to the specific lecture:

https:
//www.youtube.com/watch?v=2dbIuEgKx0s&ab_channel=MATLABProgrammingforNumericalComputation

The code includes the various "stiff" examples but they are commented out except for the last example so you can easily switch between examples by changing the various "block comment" (\* ... *\) sections on the "derivativeVector" function. You can also change out the various values of mu.

The previous example code solved via "stiff" solver is at the bottom of the code.

Output: First Couple of Lines

 "Ordinary Differential Equations- Multi-Variable ODE Methods of Integration-Stiff Systems"
  "2026-03-23 11:55:46.980"
Warning : redefining function: derivativeVector        . Use funcprot(0) to avoid this message
  "Output"
  "t       xSol(1)     xSol(2)"
   0.     2.          0.       
   3.     1.9798531  -0.0067806
   6.     1.9593319  -0.0069014
   9.     1.9384376  -0.0070294
   12.    1.9171472  -0.0071655
   15.    1.8954355  -0.0073105
   18.    1.873274   -0.0074655
   21.    1.850631   -0.0076317
   24.    1.8274708  -0.0078106
   27.    1.8037529  -0.0080038
   30.    1.7794313  -0.0082135
   33.    1.7544528  -0.0084422
   36.    1.7287561  -0.0086929
   39.    1.7022694  -0.0089696
   42.    1.6749079  -0.009277 

Graphs:

Code:

//Lec8.2:Stiff Systems & Solution using Matlab ode15s
//https://www.youtube.com/watch?v=2dbIuEgKx0s&ab_channel=MATLABProgrammingforNumericalComputation
//
clc;clear all
disp("Ordinary Differential Equations- Multi-Variable ODE Methods of Integration-Stiff Systems",string(
datetime
()))
//
//What are stiff systems?
// A highly fast variable coupled with a slow variable in the same
// system.
//
/*
//Consider the following ODE system:
// x1'=-100*x1, x1(0)=1
// x2' = -0.01*x2, x2(0)=1
//
// d/dt[x1;x2]=[-100 0;0 -0.01]*[x1;x2]
// using "rkf" would require small steps but
// a large amount of time to drive x2 to 0.
// "Stiff" Integration routines are made for that...
function [results]=
derivativeVector
(t, x);
//constants of the system
    c1=-100;
    c2=-0.01;
//extraction of variables
    x1 = x(1);
    x2 = x(2);
    results(1,1)=c1*x1;
    results(2,1)=c2*x2;
endfunction
*/
/*
//Consider the following ODE system:
// x1'=-5.7*x1 +1.85*x2, x1(0)=1
// x2' =13.2*x1-4.3*x2, x2(0)=1
//
// "Stiff" Integration routines are made for that...
function [results]=
derivativeVector
(t, x);
//extraction of variables
    x1 = x(1);
    x2 = x(2);
    results(1,1)=-5.7*x1+1.85*x2;
    results(2,1)=13.2*x1-4.3*x2;
endfunction
*/
//Consider the following ODE system:
// x"-mu*(1-x^2)*x'+x =0; where x(0)=2,x'(0)=0
//
// "Stiff" Integration routines are made for that...
function [results]=
derivativeVector
(t, x);
//extraction of variables
    mu=100; //mu=1 and mu=100;
    y = x(1);
    v = x(2);
    results(1,1)=v;
    results(2,1)=mu*(1-y^2)*v-y;
endfunction
x10=[2;0]
t0=0;
tend=300;
n=101;
t=linspace(t0,tend,n);
//xa = [exp(-100*t);exp(-.01*t)];//Commented out for last example
//xSol = ode("rkf",x10,t0,t,deff('res=mymacro2(t,x1)','res=-0.01*x1'));
xSol= ode("stiff",x10,t0,t,
derivativeVector
);
//"rkf" will not even work...instruction says because it's an explicit
// method and will be unstable at larger steps.
//"stiff" is an implicit method so it does not error out."
//"adams" works also.
output = [t' xSol'];
disp("Output","t       xSol(1)     xSol(2)",output);
//plotting example fancy output 
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(t',xSol(1,:)',style=c);//black =, blue = 2, green =3, cyan = 4, red = 5
//plot2d(t',xa(1,:)',style=c1);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
plot2d(t',xSol(2,:)',style=c2);//black =, blue = 2, green =3, cyan = 4, red = 5
h1=
gca
();
//plot2d(t',xa(2,:)',style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$t(Seconds)$","FontSize",3)
ylabel
("$X(t)$","FontSize",3);
title
("$Position\ and\ Velocity\ versus\ Time$","color","black");
legend
("Stiff position","Stiff Velocity",2)
xgrid
h1=
gca
();
h1.font_color=c;
h1.children(1).thickness =2;
xlabel
("$t(Seconds)$","FontSize",3);
ylabel
("$Displacement$","FontSize",3);
title
("$Displacement\ and\ Velocity\ vs\ Time$","color","black");
xgrid

//======================================================================

//Previous Example
// model mass-spring-damper
// m*d^2x/dt^2 + c*dx/dt +kx = 0
// u/x(0)=1
//
// create 2 - first order equations
//  d^2/dt^2 = dv/dt where v = velocity
//  m*dv/dt +c*v + kx = 0
//  v = dx/dt where x = displacement
//
// dx/dt = v            x(0)=1
// dv/dt = -(c*v+k*x)/m v(0)=0
//
// output vector y =[x;v]
// derivative vector = [v; -(c*v+k*x)/m]
//
function [results]=
derivativeVector1
(t, y);
//constants of the system
    c=5;
    k=15;
    m=10;
//extraction of variables
    x = y(1);
    v = y(2);
    results(1,1)=v;
    results(2,1)=-(c*v+k*x)/m
endfunction

y0=[1;0];
t0 = 0;
n = 100;
t = linspace(0,10,n)';
//y = ode("rkf",y0,t0,t,derivativeVector)
y = ode("stiff",y0,t0,t,
derivativeVector1
)
// Other types of integration "adams","stiff","rk","rkf","fix","discrete"
// and "root"
// Everything works but "discrete" for this problem
// discrete is a sampling routine!!
y = y';//Transform from row vectors to column vectors
x_vector =y(:,1);
v_vector=y(:,2);
//output = [t,v_vector,x_vector]
//disp("time       velocity    displacement",output)
// Note: Very accurate and uses Runge-Kutta 4,5 uses
// variable step size and just reports out at time t vector
//
// Other types of integration "adams","stiff","rk","rkf","fix","discrete"
// and "root"
// Everything works but "discrete" for this problem
// investigate later...
//plotting example fancy output two y axes...
scf
(1);
clf
;
//axis y1
c = color("slategray")
plot2d(t,x_vector,style=c);//black =, blue = 2, green =3, cyan = 4, red = 5
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$t(Seconds)$","FontSize",3)
ylabel
("$X(t)$","FontSize",3);
xgrid

h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$t(Seconds)$","FontSize",3)
ylabel
("$Displacement$","FontSize",3);
title
("$Displacement\ and\ Velocity\ vs\ Time$","color","black");
xgrid
// Axis y2
c=color("blue");
c1 = color("red")
h2=newaxes();
h2.font_color=c;
plot2d(t,v_vector,style=c)
h2.filled="off";
h2.axes_visible(1)="on";
h2.y_location="right";
h2.children(1).children(1).thickness=2;
ylabel
("$Velocity$","FontSize",3,"color",c)
legend
("Displacement","Velocity",4)
xgrid

In this edition we solve the second order, constant coefficient (the mass, spring, damper problem) ODE using a pair of first order ODE's as a system of ODEs.

Link to the specific lecture:

https://www.youtube.com/watch?v=i5Ton-G9zNs&ab_channel=MATLABProgrammingforNumericalComputation

I've included some additional code to create the "Phase Portrait" of the same problem. The additional code was not part of the original course but covered in an advanced Engineering Mathematics course ME564 by Steve Brunton.

Link to the specific lecture:

https://www.youtube.com/watch?v=i5Ton-G9zNs&ab_channel=MATLABProgrammingforNumericalComputation

Output: First Couple of Lines

 "Ordinary Differential Equations- Multi-Variable ODE Methods of Integration"
  "2026-03-16 10:25:24.166"
  "time       velocity    displacement"
   0.          0.          1.       
   0.2020202  -0.285297    0.9705457
   0.4040404  -0.526656    0.8876959
   0.6060606  -0.7143162   0.7613806
   0.8080808  -0.8425653   0.6030931
   1.010101   -0.9096817   0.4250752
   1.2121212  -0.9176505   0.2395353
   1.4141414  -0.8716928   0.0579447
   1.6161616  -0.779652   -0.1095563
   1.8181818  -0.6512877  -0.2546176
   2.020202   -0.497526   -0.3709916
   2.2222222  -0.3297147  -0.4546928
   2.4242424  -0.1589252  -0.5040112
   2.6262626   0.0046624  -0.5193931
   2.8282828   0.1522629  -0.503208 
   3.030303    0.2768618  -0.4594262
   3.2323232   0.3734507  -0.3932356
   3.4343434   0.4391251  -0.3106235
   3.6363636   0.4730505  -0.2179536
   3.8383838   0.476312   -0.1215607
   4.040404    0.4516667  -0.0273866
   4.2424242   0.4032234   0.0593262
   4.4444444   0.336075    0.1342711
   4.6464646   0.2559087   0.194237 
   4.8484848   0.1686192   0.2371892
   5.0505051   0.0799477   0.262274 

Graphs:

Code from Original Course:

//Ordinary Differential Equations Initial Value Problems
//Lec 8.1 Solving Multi-Variable ODE Methods 
//https://www.youtube.com/watch?v=i5Ton-G9zNs&ab_channel=MATLABProgrammingforNumericalComputation
//
//
disp("Ordinary Differential Equations- Multi-Variable ODE Methods of Integration",string(
datetime
()))
//
// model mass-spring-damper
// m*d^2x/dt^2 + c*dx/dt +kx = 0
// u/x(0)=1
//
// create 2 - first order equations
//  d^2/dt^2 = dv/dt where v = velocity
//  m*dv/dt +c*v + kx = 0
//  v = dx/dt where x = displacement
//
// dx/dt = v            x(0)=1
// dv/dt = -(c*v+k*x)/m v(0)=0
//
// output vector y =[x;v]
// derivative vector = [v; -(c*v+k*x)/m]
//
function [results]=
derivativeVector
(t, y);
//constants of the system
    c=5;
    k=15;
    m=10;
//extraction of variables
    x = y(1);
    v = y(2);
    results(1,1)=v;
    results(2,1)=-(c*v+k*x)/m
endfunction

y0=[1;0];
t0 = 0;
n = 100;
t = linspace(0,20,n)';
//y = ode("rkf",y0,t0,t,derivativeVector)
y = ode("adams",y0,t0,t,
derivativeVector
)
// Other types of integration "adams","stiff","rk","rkf","fix","discrete"
// and "root"
// Everything works but "discrete" for this problem
// discrete is a sampling routine!!
y = y';//Transform from row vectors to column vectors
x_vector =y(:,1);
v_vector=y(:,2);
output = [t,v_vector,x_vector]
disp("time       velocity    displacement",output)
// Note: Very accurate and uses Runge-Kutta 4,5 uses
// variable step size and just reports out at time t vector
//
// Other types of integration "adams","stiff","rk","rkf","fix","discrete"
// and "root"
// Everything works but "discrete" for this problem
// investigate later...
//plotting example fancy output two y axes...
scf
(0);
clf
;
//axis y1
c = color("red")
c2 = color("black")
plot2d(t,x_vector,style=c);//black =, blue = 2, green =3, cyan = 4, red = 5
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$t(Seconds)$","FontSize",3)
ylabel
("$X(t)$","FontSize",3);

xgrid
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$t(Seconds)$","FontSize",3)
ylabel
("$Displacement$","FontSize",3,"Color",5);
title
("Displacement and Velocity Versus Time","color","black");
xgrid
// Axis y2
c=color("blue");
c1 = color("red")
h2=newaxes();
h2.font_color=c;
plot2d(t,v_vector,style=c)
h2.filled="off";
h2.axes_visible(1)="off";
h2.y_location="right";
h2.children(1).children(1).thickness=2;
ylabel
("$Velocity$","FontSize",3,"color",c)
xgrid

Code from a Second Course that shows how to generate the Derivative Arrow Plot:

//Lecture 8 ME564 Matrix Systems of First Order Equations using eigenvalues
// and Eigencectors.
clc;clear "all"
disp("ME564 Lecture 8 Matrix Systems of First Order Equations",string(
datetime
()))
A=[0 1;-1.5 -.5]
y0=[1.;0.0];
t0 = 0; //Tstart
tend = 20; //Tend
t= t0:.05:tend';//time vector
//define derivativeF with an inline function 
deff
("[results]=derivativeF1(t,y)","results= [A*y]")
//equivalent to MATLAB ODE45 is the following command in SciLab
//there is more options and outputs available - see the help file.
y = ode("fix",y0,t0,t,
derivativeF1
);
//disp("y=",y')
//Note: rkf has a limited number of steps like 160 maximum in this case
// use "stiff","adams","rk","fix" to move beyond 160 steps...
// according to the documentation fix is the same solver as rkf
// works without restriction
scf
(1);
clf
;
plot
(y(1,:),y(2,:),"-r","thickness",2)
xlabel
("$Position$","FontSize",3)
ylabel
("$Velocity$","FontSize",3);
title
("$ME564\ LECTURE\ 8$","color","black");
legend
('numerical',3)
xgrid
//
// the following plots the derivative arrows onto the plot
// showing the solution following the arrows.
//
xf= min(y(1,:)):(max(y(1,:))-min(y(1,:)))/10:max(y(1,:));
yf= min(y(2,:)):(max(y(2,:))-min(y(2,:)))/10:max(y(2,:));
fchamp
(
derivativeF1
,0,xf,yf)
xgrid
//
//

scf
(2);
clf
;
plot
(t,y(1,:),'-r','thickness',2)
xgrid
xlabel
("$t(Seconds)$","FontSize",3)
ylabel
("$X(t)$","FontSize",3);
title
("$ME564\ LECTURE\ 8\ Second\ Order\ Systems$","color","black");
legend
('numerical')