Euler Method System of ODE solving

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);

Answer 1

[Alan Stevens]

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')

[Ajmal R S]

Thank you, the code is working now, though there seems to be some error n function 2, temperature is not supposed to hit values of 10^303, I might have made some error with that function

[Ajmal R S]

This was the second function in the euler method code that I did, but when i plot the graph I am getting values of T2 to around 4*10^303, which is wrong. I tried solving the equation using Ode15s and got the correct answer of 318.5 K, please help

[Alan Stevens]

Simple Euler is well known to be inaccurate, especially for complicated equations. Never use it for anything other than "toy" problems. MATLAB's built-in solvers are the way to go. Since you used ode15s to get a good solution I assume it's a stiff problem, which makes it even more unlikely that simple Euler will get a sensible result.

(Incidentaly, you didn't say what your initial conditions or final time are).

[Ajmal R S]

Initial condition for T2 is 320.1K and the time span is from [0,200]. I tried to use Simple euler to compare the steady state value with the steady state value obtained from the ode15s and to yield the SSE values