Solving system of ODEs with ode45

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

[Walter Roberson]

What is the evidence that those are not correct output graphs?

Answer 1

[Sam Chak]

Hi @shir levi

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

Answer 2

[shir levi]

thak you for answaring

but those arnt the initial conditions

its the valuse i was supposed to get in my results

[Sam Chak]

Hi @shir levi

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.