# Solving system of ODEs with ode45

1 view (last 30 days)
shir levi on 24 Dec 2022
Commented: Sam Chak on 25 Dec 2022
i have the following equetions system
and and the maching valuse
i tried to solve the system with the following code but i dont get the right resultes which i know:
ephsilon values [1.94 2.65 3.42] for the x1/h vec [7.5 9.5 11]
k values [0.2680 0.3585 0.444] for the same x1/h vec.
the code i wrote
Y=zeros(1,2);
Uc=12.4; %m/s
S=46.8;% 1/s
k_75=0.268 ; %m2/s2;
e_75=1.94; %m2/s3
Cu = 0.09;
C1= 1.44;
C2=1.92;
x1_h=7.5;
h = 0.3048;%m
xspan=[7.5 9.5 11];
Y0=[k_75 e_75 ];
f = @(t,Y) [(Cu*((Y(1))^2)/(Y(2))*(S^2)-Y(2))/(Uc/h); (Cu*C1*(Y(1))*(S^2)-C2*((Y(2))^2)/(Y(1)))/(Uc/h)];
[x0,val] = ode45(f , xspan, Y0);
figure (1)
plot(x0, val(:,1))
figure (2)
plot(x0, val(:,2))
what wrong with it ? thanks for the help
##### 1 CommentShowHide None
Walter Roberson on 25 Dec 2022
What is the evidence that those are not correct output graphs?

Sam Chak on 24 Dec 2022
Edited: Sam Chak on 24 Dec 2022
Edit: Previously I've got the same results like you did. However, after thinking about it, those 2 sets of 3 values for k and ϵ are probably initial conditions. So, I modified the code to include a for loop.
Also, in the 3rd figure, the system is shown to be inherently unstable, and there is no equilibrium point.
% parameters
Uc = 12.4; % m/s
S = 46.8; % 1/s
k = [0.2680 0.3585 0.444]; % m2/s2 ic for k
epsil = [1.9400 2.6500 3.420]; % m2/s3 ic for epsil
Cu = 0.09;
C1 = 1.44;
C2 = 1.92;
x1_h = [7.5 9.5 11];
h = 0.3048; % m
tstart = x1_h*h/Uc;
tfinal = 2*tstart;
for j = 1:length(x1_h)
% setting conditions
f = @(t, Y) [Cu*(Y(1)^2)/Y(2)*S^2 - Y(2); Cu*C1*Y(1)*S^2 - C2*(Y(2)^2)/Y(1)];
tspan = [tstart(j) tfinal(j)];
Y0 = [k(j) epsil(j)];
[t, Y] = ode45(f, tspan, Y0);
% plotting results
figure(1)
plot(t*Uc/h, Y(:,1)), hold on, grid on,
xlabel('x_{1}/h'), ylabel('k'), title('k responses')
figure(2)
plot(t*Uc/h, Y(:,2)), hold on, grid on,
xlabel('x_{1}/h'), ylabel('\epsilon'), title('\epsilon responses')
end
hold off
[x, y] = meshgrid(-11:2:11);
u = Cu*(x.^2)./y*S^2 - y;
v = Cu*C1*x*S^2 - C2*(y.^2)./x;
figure(3)
l = streamslice(x, y, u, v);
axis tight
xlabel('k')
ylabel('\epsilon')

shir levi on 24 Dec 2022
thak you for answaring
but those arnt the initial conditions
its the valuse i was supposed to get in my results
##### 1 CommentShowHide None
Sam Chak on 25 Dec 2022
But you used these values and in the initial condition.
So, I'm confused.
Perhaps, if you tell us more info about the 2 sets of 3 values, it will help clarifying the matter.