I am having trouble multiplying my main ode function with an external function NS which is supposed to be multiplied on the RHS. Thanks for the great help.

tspan = [0 20];
y0 = [0 0.01];
[z,y] = ode45(@odefcn, tspan, y0);
plot(z,y(:,1),'-o',z,y(:,2),'-.')

function dydt = odefcn(z,y)
dydt1 = y(2);
dydt2 = NS*z.*y(1);
dydt=[dydt1;dydt2];
end 

function M = NS(z)

z = [2 3 5 7 10 15 20 ];
r =[3.5 3.7 4 6 7.2 8 9];
n=length(z);
% Calculation of differentiation from the above datas
  for i=1:n-1
     M=zeros();
    M(i)=(r(i+1)-r(i))./(z(i+1)-z(i));

    end
  end

[Torsten]

It's not obvious which ODEs you are trying to solve.

What has "t" to do with "r" and "z" ?

Best wishes

Torsten.

[Dereje]

yess you are right Torsten, now I corrected to z(I was supposed to mean z not t)

[Dereje]

I wanted to solve a function with the name odefcn. But my problem is the RHS is supposed to be multiplied with external function which is differentiation of data points.

Answer 1

[Torsten]

zr = [2 3 5 7 10 15 20];
r = [3.5 3.7 4 6 7.2 8 9];
y0 = [0 ; 0.01];
zsol = [];
y1sol = [];
y2sol = [];
for i=1:numel(r)-1
  zspan = [zr(i) zr(i+1)];
  NS = (r(i+1)-r(i))/(zr(i+1)-zr(i));
  [z,y] = ode45(@(z,y)odefcn(z,y,NS), zspan, y0);
  y0 = [y(end,1) ; y(end,2)];
  zsol = [zsol;z];
  y1sol = [y1sol;y(:,1)];
  y2sol = [y2sol;y(:,2)];
end
plot(zsol,y1sol,zsol,y2sol)    

function dydt = odefcn(z,y,NS)
dydt1 = y(2);
dydt2 = NS*z*y(1);
dydt=[dydt1;dydt2];
end

Best wishes

Torsten.

[Dereje]

Yes!! It worked.Thank you very much!

[Dereje]

How can I correct for the pre-allocation error that highlights in red on :

 zsol =  [zsol;t];
v1sol = [v1sol;v(:,1)];
v2sol = [v2sol;v(:,2)];
v3sol = [v3sol;v(:,3)];

[Torsten]

Please show your code and the error message you receive.

Best wishes

Torsten.

[Dereje]

I have attached the code and the highlight is on

 zsol =  [zsol;t];
v1sol = [v1sol;v(:,1)];
v2sol = [v2sol;v(:,2)];
v3sol = [v3sol;v(:,3)];

Thanks a lot!

[Torsten]

I must admit that I don't see an error in the code.

Does this work for you ?

%%Data for desnity with respect to depth 

z = [2 3 5 7 10 15 20 25 30 40 50 60 70 80 90 100 125 150 160 175 200 225 250 275 300 325 350 375 400];
rho = [17.2731684 17.1649375 21.43455647 22.4140625 23.86332207 24.3746967 24.70487685 24.6003125 24.8933125 25.42772826 26.03220776 26.439625 26.8151875 26.86830797 27.1949375 27.34406944 27.5551875 27.728625 27.23423729 27.88542857 27.752249049 28.1025 28.2415 28.37 28.05366667 28.6565 28.7755 28.898 29.013]; 

v0=[0.8;0.001;0.1;]; % Initial values 

% Creating a matrix 
zsol=zeros(10000);
v1sol=zeros(10000);
v2sol=zeros(10000);
v3sol=zeros(10000);
start=1;
for i=1:numel(rho) - 1
rho0=17.1;
g=9.8;
zspan = [z(i) z(i+1)];
Nsquare = (g/rho0)*(rho(i+1)-rho(i))/(z(i+1)-z(i));
[t,v] = ode45(@(t,v)rhs(t,v,Nsquare), zspan, v0);
ende = start+numel(t);
v0 = [v(end,1) ; v(end,2) ; v(end,3)];
zsol(start:ende)=t(:);
v1sol(start:ende)=v(:,1);
v2sol(start:ende)=v(:,2);
v3sol(start:ende)=v(:,3);
start=ende;
end
plot(zsol(1:ende),v1sol(1:ende),'r',zsol(1:ende),v2sol(1:ende),'g',zsol(1:ende),v3sol(1:ende),'b')
xlabel('Width,b and vertical velocity,w')
ylabel('Height, z')
grid on ;
function parameters=rhs(t,v,Nsquare)

       alpha=0.116;
              db= 4*alpha*v(2).^2-v(1).*v(3)./2*v(2).^2;
              dw= v(1).*v(3)-4*alpha*v(2).^4+v(1).*v(2).^2.*v(3)./2*v(1).*v(2).^3;
              dgmark= (-2*Nsquare*v(1).*v(2)^4-v(1).*v(3)^2+4*alpha*v(2).^4.*v(3)-v(1).*v(2).^2.*v(3).^2-8*alpha*v(2).^3.*v(3)+2*v(3).^2.*v(1).*v(2))./2*v(1).*v(2).^4; 

     parameters=[db;dw;dgmark];
   end

[Dereje]

Now the highlight(pre-allocation) is gone but when I run it(the one you corrected above) It has an error

(Error line 20)
zsol(start:ende)=t(:);

In an assignment

 A(:) = B, the
number of elements in A and B must
be the same.

[Torsten]

Use

ende = start+numel(t)-1;

instead of

ende = start+numel(t);

Best wishes

Torsten.

[Dereje]

It works but it takes longer time to execute the code. one time it crushed. Now I am trying to avoid this. Do you have any tips ? Thanks

[Torsten]

Using ode15s instead of ode45 might reduce the computation time.

Best wishes

Torsten.

[Dereje]

ya a little better! Thanks!!