Solving the differential equation by the Runge-Kutta method (ode45)

I want solve this equation numerically and use fourth order Runge-Kutta method. (ode45)

What is the code?

My code:

function dfdt = odefun2(x,y)
dfdt = [y(2); y(1) .^(-3) - (w .*r0 ./c) .^2 .*phi2 .*y(1)];
end

****

clc, clear, close all;
w = 2e16;
r0 = 30e-6;
c = 3e8;
p = 2.5;
f = 1;
Wp0 = p .*c ./r0;
a02 = 0.050;
phi2 = (Wp0 ./w) .^2 .*(3 .*a02 ./f .^4) ./4 ./(1 + a02 ./2 ./f .^2) .^1.5...
    .*(1 + (56 + a02 ./f .^2) ./3 ./(1 + a02 ./2 ./f .^2) ./p .^2 ./f .^2);

[t, y] = ode45(@odefun2,[0 3],[1,0]);

******

But it gives an error ..

I have attached two pictures to the question, look at them if you need.

[Aquatris]

The expectation on this site is for you to try it yourself first and ask questions when you get stuck. So what did you try and where you got stuck?

[Torsten]

What boundary conditions for f ?

[Babr]

w = 2e16;

r0 = 30e-6;

c = 3e8;

p = 2.5;

f = 1;

Wp0 = p .*c ./r0;

a02 = 0.050;

phi2 = (Wp0 ./w) .^2 .*(3 .*a02 ./f .^4) ./4 ./(1 + a02 ./2 ./f .^2) .^1.5...

.*(1 + (56 + a02 ./f .^2) ./3 ./(1 + a02 ./2 ./f .^2) ./p .^2 ./f .^2);

[Torsten]

According to your code, your boundary conditions are

f(zeta=0) = 1
df/dzeta (zeta=0) = 0

Is this correct ?

Did you differentiate depsilon/dzeta somewhere to insert it into the differential equation ? I cannot find it in your code.

[Babr]

Yes, the boundary conditions you said are correct.

The second term in right hand of equation has been usually overlooked in most of the studies in view of negligible impact.

Answer 1

[Torsten]

[t, y] = ode45(@odefun2,[0 3],[1,0]);
plot(t,y(:,1))

function dfdt = odefun2(x,y)
  w = 2e16;
  r0 = 30e-6;
  c = 3e8;
  p = 2.5;
  Wp0 = p .*c ./r0;
  a02 = 0.050;
  phi2 = (Wp0/w)^2 * 3*a02/y(1)^4 / (4*(1+a02/(2*y(1)^2))^1.5) *...
         (1 + 1/(p^2*y(1)^2) * (56+a02/y(1)^2)/(3*(1+a02/(2*y(1)^2))^0.5));   
  dfdt = [y(2); y(1) .^(-3) - (w .*r0 ./c) .^2 .*phi2 .*y(1)];
end

[Babr]

Thank you so much.

I have another question.

If we want to consider the second term on the right hand of the equation, then how can the equation be solved using the ode command?

[Torsten]

Take paper and pencil and differentiate the expression for epsilon0' with respect to zeta. Then insert the result into the differential equation as you did for phi2.

Or do it symbolically with MATLAB and copy the resulting expression into your code.

[Babr]

Thank you very much for your guidance and help and for informing me of my mistake.

I wish you success and happiness.

[Torsten]

syms f(x) p c r0 a02 w
Wp0 = p*c/r0;
e = 1 - (Wp0/w)^2 / sqrt(1+a02/(2*f^2)) * (1- 1/p^2 * 3*a02/f^4 / sqrt(1+a02/(2*f^2)));
simplify(diff(e,x))

ans(x) = 

and df/dx in your code is y(2).

[t, y] = ode45(@odefun2,[0 3],[1,0]);
plot(t,y(:,1))

function dfdt = odefun2(x,y)
  w = 2e16;
  r0 = 30e-6;
  c = 3e8;
  p = 2.5;
  Wp0 = p .*c ./r0;
  a02 = 0.050;
  phi2 = (Wp0/w)^2 * 3*a02/y(1)^4 / (4*(1+a02/(2*y(1)^2))^1.5) *...
         (1 + 1/(p^2*y(1)^2) * (56+a02/y(1)^2)/(3*(1+a02/(2*y(1)^2))^0.5));
  e = 1 - (Wp0/w)^2 / sqrt(1+a02/(2*y(1)^2)) * (1- 1/p^2 * 3*a02/y(1)^4 / sqrt(1+a02/(2*y(1)^2)));
  sigma1 = a02/(2*y(1)^2) + 1;
  dedx = 3*a02^2*c^2*y(2)/(r0^2*w^2*y(1)^7*sigma1^2) ...
         -12*a02*c^2*y(2)/(r0^2*w^2*y(1)^5*sigma1)...
         -a02*c^2*p^2*y(2)/(2*r0^2*w^2*y(1)^3*sigma1^1.5);
  dfdt = [y(2); y(1)^(-3)-1/(2*e)*dedx*y(2)-(w*r0/c)^2*phi2*y(1)];
end

[Babr]

I don't know how to thank you ..