MATLAB: Euler Method System of ODE solving

MATLAB and Simulink Student Suiteode

I have tried to solve a system of 3 ODEs , and this is my code
clear all
clc
V1  = 1.75e-03;   %Volume Tank 1
Tc  = 303;        %Cooling Water Temp

T2_a  = 320.1;      %Tank 2 Temp
rho_w = 1000;     %Density in kg/m3

Cp  = 4182;       %Cp value in J/KgK

u1  = 45;       %Flow F1


F1_A =  8.0368e-11;
F1_B =  456.85e-11;
F1_C =  42379e-11;
F1 = F1_A*(u1^3)-F1_B*(u1^2)+F1_C*(u1);
u3  = 45;       %Flow F3

Fr_A =  (1/1800)*1e-03;
Fr = Fr_A*u3;
u4  = 55;       %Heat Input Q1
Q1_A =  7.9798;
Q1_B =  0.9893;
Q1_C =  0.0073;
Q1 = Q1_A*(u4)+Q1_B*(u4^2)-Q1_C*(u4^3);
%dy(1) = (F1*(Tc-y(1))+Fr*(T2_a-y(1))+(Q1/(rho_w*Cp)))/(V1);
A2  = 7.854e-05;  %Cross-Sectional Area Tank 2

Tc  = 303;        %Cooling Water Temp
T1_a  = 325.4;      %Tank 1 Temp
Ta  = 298;        %Atmospheric Temp
h2  = 0.41;       %Level h2
R2  = 0.05;       %Radius Tank 2
rho_w = 1000;     %Density in kg/m3
Cp  = 4182;       %Cp value in J/KgK
U   =  235.1043;   %Heat Transfer Coeff
u1  = 45;       %Flow F1
F1_A =  8.0368e-11;
F1_B =  456.85e-11;
F1_C =  42379e-11;
F1 = F1_A*(u1^3)-F1_B*(u1^2)+F1_C*(u1);
u2  = 45;       %Flow F2

F2_A =  1.2943e-11;
F2_B =  190.64e-11;
F2_C =  8796.8e-11;
F2_D =  196620e-11;
F2 = -F2_A*(u2^4)+F2_B*(u2^3)-F2_C*(u2^2)+F2_D*(u2);
u3  = 45;       %Flow F3
Fr_A =  (1/1800)*1e-03;
Fr = Fr_A*u3;
u5  = 60;       %Heat Input Q2
Q2_A = 14.4;
Q2_B = 0.96;
Q2_C = 0.008;
Q2 = Q2_A*(u5)+Q2_B*(u5^2)-Q2_C*(u5^3);
%dy(2) = (F1*(T1_a-y(2))+F2*(Tc-y(2))-Fr*(y(2)-T1_a)+(Q2-(2*pi*h2*R2*U*(y(2)-Ta))/rho_w*Cp))/(A2*h2);     
A2  = 7.854e-05;  %Cross-Sectional Area Tank 2
k   =  0.1;        %Discharge Coefficient
u1  = 45;       %Flow F1
F1_A =  8.0368e-11;
F1_B =  456.85e-11;
F1_C =  42379e-11;
F1 = F1_A*(u1^3)-F1_B*(u1^2)+F1_C*(u1);
u2  = 45;       %Flow F2
F2_A =  1.2943e-11;
F2_B =  190.64e-11;
F2_C =  8796.8e-11;
F2_D =  196620e-11;
F2 = -F2_A*(u2^4)+F2_B*(u2^3)-F2_C*(u2^2)+F2_D*(u2);
Fd_A = 0.4060;
Fd_B = 0.80608;
Fd_C = -0.01798;
Fd_D = 0.10541;
%dy(3) = (F1+F2-k*((Fd_A*((y(3)^3))+Fd_B*((y(3)^2))-Fd_C*(y(3))+0.10541)^0.5)*1e-03)/A2;
f=@(t,y)([([(F1*(Tc-y(1))+Fr*(T2_a-y(1))+(Q1/(rho_w*Cp)))/(V1);(F1*(T1_a-y(2))+F2*(Tc-y(2))-Fr*(y(2)-T1_a)+(Q2-(2*pi*h2*R2*U*(y(2)-Ta))/rho_w*Cp))/(A2*h2); (F1+F2-k*((Fd_A*((y(3)^3))+Fd_B*((y(3)^2))-Fd_C*(y(3))+0.10541)^0.5)*1e-03)/A2;])]);    
 %Write your f(x,y) function, where dy/dx=f(x,y), x(x0)=y0.
x0=input('\n Enter initial value of x i.e. x0: '); %example x0=0
y0=input('\n Enter initial value of y i.e. y0: '); %example y0=0.5
xn=input('\n Enter the final value of x: ');% where we need to find the value of y 
                                            %example x=2
h=input('\n Enter the step length h: '); %example h=0.2
 %Formula: y1=y0+h*fun(x0,y0);
 fprintf('\n x        y ');
 
 
while x0<=xn
    i=1;
    fprintf('\n%4.3f  %4.3f ',x0,y0); %values of x and y
    y1(i)=y0+h*f(x0,y0);
    x1=x0+h;
    x0=x1;
    y0=y1(i);  
    i=i+1;
end
But I got this error, please advice
Index exceeds the number of array elements (1).
Error in
euler_csth>@(t,y)([([(F1*(Tc-y(1))+Fr*(T2_a-y(1))+(Q1/(rho_w*Cp)))/(V1);(F1*(T1_a-y(2))+F2*(Tc-y(2))-Fr*(y(2)-T1_a)+(Q2-(2*pi*h2*R2*U*(y(2)-Ta))/rho_w*Cp))/(A2*h2);(F1+F2-k*((Fd_A*((y(3)^3))+Fd_B*((y(3)^2))-Fd_C*(y(3))+0.10541)^0.5)*1e-03)/A2])])
(line 99)
f=@(t,y)([([(F1*(Tc-y(1))+Fr*(T2_a-y(1))+(Q1/(rho_w*Cp)))/(V1);(F1*(T1_a-y(2))+F2*(Tc-y(2))-Fr*(y(2)-T1_a)+(Q2-(2*pi*h2*R2*U*(y(2)-Ta))/rho_w*Cp))/(A2*h2);
(F1+F2-k*((Fd_A*((y(3)^3))+Fd_B*((y(3)^2))-Fd_C*(y(3))+0.10541)^0.5)*1e-03)/A2;])]);
Error in euler_csth (line 113)
y1(i)=y0+h*f(x0,y0);

Best Answer

Your function, f = @(t,y) ... has vector y with just three components, but you call it in your Euler routine with i = 4 ...
y needs to be a matrix of 3xn where n is the number of points at which you want to evauate the function (which itself had a syntax error).
If you replace all the lines in your code, from that function definition, with the following code, your program will work. (You will, of course, have to replace the arbitrary constants/times/etc that I've used with your actual data).
f=@(t,y)[(F1*(Tc-y(1))+Fr*(T2_a-y(1))+(Q1/(rho_w*Cp)))/(V1);
       (F1*(T1_a-y(2))+F2*(Tc-y(2))-Fr*(y(2)-T1_a)+(Q2-(2*pi*h2*R2*U*(y(2)-Ta))/rho_w*Cp))/(A2*h2);
       (F1+F2-k*((Fd_A*((y(3)^3))+Fd_B*((y(3)^2))-Fd_C*(y(3))+0.10541)^0.5)*1e-03)/A2];    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following data is arbitrary and needs to be replaced by true values
t0 = 0;  tf = 1; % Initial and final times/positions respectively
n = 100;         % Number of times/segments
h = (tf-t0)/n;   % step size
t = t0:h:tf;     % times
y = zeros(3,n+1); % Storage space
y(1,1) = Tc; y(2,1) = T1_a; y(3,1) = 1;  % Initial values
% Simple Euler integration
for i = 1:n
    
    y(:,i+1)=y(:,i)+h*f(t(i),y(:,i));
    
end
% Plot results 
subplot(3,1,1)
plot(t,y(1,:)),grid
xlabel('t'),ylabel('y1')
subplot(3,1,2)
plot(t,y(2,:)),grid
xlabel('t'),ylabel('y2')
subplot(3,1,3)
plot(t,y(3,:)),grid
xlabel('t'),ylabel('y3')
Related Question