solving the equation using rk-4 method, but couldn't get the right output, last time I forget to add the program please check this one.

i want to solve this, where σ = 0.1, L0 = 0.5, k= 0.01, tc = 70E-9

format long

a = 0;

b = 4;

h = 0.1;

k = 0.01; % critical coupling strength

L = 0.5;

s = 0.1;

tc = 70E-9; % photon life time in cavity

n = (b-a)/h ;

t = a + (0:n)*h;

func1 = @(t,y) [(-(s^2)*k/tc)*sin(y - pi/2) + L*(s^2)/(2*tc)*sin(y + pi/2)/sqrt(1 + cos(y + pi/2))];

for i = 1:n

k1 = func1(y(i));

k2 = feval(func1, y(i)+(k1/2)*h);

k3 = feval(func1, y(i)+(k2/2)*h);

k4 = feval(func1, y(i)+(k3*h));

y(i+1) = y(i) + (h/6)*(k1+(2*(k2+k3))+k4);

end

plot(t,y)

[Torsten]

You posted an algebraic equation connecting U and phi0. Why do you solve a differential equation ?

[SAHIL SAHOO]

sorry by mistake I added different equation, please check now.

[Torsten]

You don't supply an initial condition for phi0 (must be written in y(1)).

Further, your function handle depends on t and y. You try to evaluate it only for one input (y).

[SAHIL SAHOO]

it solved, thanks man

[SAHIL SAHOO]

@Torsten i made it properly still it's not working.

================================================

format long
a = 0;
b = 10E-5;
h = 0.01E-5;
k = 0.01;  % critical coupling strength
L  = 0.5;
s = 0.1;
tc  = 70E-9;  % photon life time in cavity

n = (b-a)/h ;
t = a + (0:n)*h;
y(1) = 0;
func1 = @(Y) [(-(s^2)*k/tc)*sin(Y - pi/2) + L*(s^2)/(2*tc)*sin(Y + pi/2)/sqrt(1 + cos(Y + pi/2))];         

for i = 1:n
    k1 = func1(y(i));
    k2 = feval(func1, y(i)+(k1/2)*h);
    k3 = feval(func1, y(i)+(k2/2)*h);
    k4 = feval(func1, y(i)+(k3*h));
    y(i+1) = y(i) + (h/6)*(k1+(2*(k2+k3))+k4);   
    
end

plot(t,y)

=========================================

output:-

Error using plot

Vectors must be the same length.

=============================================

i don't find any problem in code, still it's not working can you please check.

[Walter Roberson]

You should preallocate y as length n+1. With your current code if you just happen to have a longer y in your workspace you would have a problem.

[Torsten]

@SAHIL SAHOO

Works for me (see above).