MATLAB: How do i solve the errors arising from this code

coupled ode equationsode45

   function my_ode()
  %parameters used inthe simulation.
  pi= 3.143; lambda_1 = 1536 ; lambda_2 = 1538; lambda_s = 1534; lambda_0 = 1535;
  lambda_c = (lambda_1 + lambda_2)/2;
  w_0 = (2*pi)/lambda_0;
  w_1 = (2*pi)/lambda_1;
  w_2 = (2*pi)/lambda_2;
  w_s = (2*pi)/lambda_s;
  w_c = (2*pi)/lambda_c;
  B_3 =2.881*10^-2;B_4 =-1.852*10^-4;P_1 = 1;P_2 =1 ;
  gamma =122.7*10^-3;
  B_2 = B_3*(w_c - w_0)^2 + (B_4/2)*(w_c - w_0);
  delta_w_p = (w_2 - w_1)/2;
  delta_w_s = w_s - w_c;
  w_sb1 = 2*delta_w_p - delta_w_s;
  %w_sb2 = 2*delta_w_p - delta_w_s;
  delta_w_sb = w_c - w_sb1;
  r = 2*gamma*sqrt(P_1*P_2);
  A_sb10 = 0;
  A_s0 = 1;
  A_10 = 0;
  A_sb20 = 0;
% function to compute derivative of
%A = [A_s, A_1, A_sb1, A_sb2];
function dAdz = odefcn(~,A)
dA_sb1dz = (1i/24)*A_sb1*(12*B_2*(delta_w_sb^2 - delta_w_p^2)... 
- 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_sb^4 - delta_w_p^4))+1i*gamma...
*P_1*(A_sb1 + conj(A_s))+ 1i*r*(conj(A_sb1) + A_s);
dA_sdz = (1i/24)*A_s(12*B_2*(delta_w_s^2 - delta_w_p^2)... 
- 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_s^4 - delta_w_p^4))+ 1i*gamma...
*P_1*(conj(A_sb1) + A_s)+ 1i*r*(conj(A_1) + A_sb2);
dA_1dz = (1i/24)*A_1(12*B_2*(delta_w_s^2 - delta_w_p^2)... 
+ 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_s^4 - delta_w_p^4))+1i*gamma...
*P_2*(conj(A_sb2) + A_1)+ 1i*r*(conj(A_s) + A_sb1);
dA_sb2dz = (1i/24)*A_sb2(12*B_2*(delta_w_sb^2 - delta_w_p^2)... 
+ 4*B_3*(delta_w_sb^3-delta_w_p^3)+ B_4*(delta_w_sb^4 - delta_w_p^4))+1i*gamma...
*P_2*(conj(A_1) + A_sb2)+ 1i*r*(conj(A_sb1) + A_s);
   dAdz = A;
   dAdz(:) = [dA_sb1dz dA_sdz dA_1dz dA_sb2dz]; 
  end
  A0 = [A_s0 A_10 A_sb10 A_sb20];
  [z,A] = ode45(@odefcn, [0 1], A0);
   A_sb1  = A(:,1); A_1 = A(:,3);
   A_s  = A(:,2);   A_sb2 = A(:,4);
  end

Best Answer

Please check whether this is what you want:
function main
  A_sb10 = 0;
  A_s0 = 1;
  A_10 = 0;
  A_sb20 = 0;
 A0 = [A_s0 A_10 A_sb10 A_sb20];
  [z,A] = ode45(@(t,y)odefcn(t,y(1),y(2),y(3),y(4)), [0 1], A0);
   A_sb1  = A(:,1); A_1 = A(:,3);
   A_s  = A(:,2);   A_sb2 = A(:,4);
   plot(z,real(A(:,4)))
  end
function dAdz = odefcn(t ,A_s, A_1, A_sb1, A_sb2)
  %parameters used inthe simulation.
  pi= 3.143; lambda_1 = 1536 ; lambda_2 = 1538; lambda_s = 1534; lambda_0 = 1535;
  lambda_c = (lambda_1 + lambda_2)/2;
  w_0 = (2*pi)/lambda_0;
  w_1 = (2*pi)/lambda_1;
  w_2 = (2*pi)/lambda_2;
  w_s = (2*pi)/lambda_s;
  w_c = (2*pi)/lambda_c;
  B_3 =2.881*10^-2;B_4 =-1.852*10^-4;P_1 = 1;P_2 =1 ;
  gamma =122.7*10^-3;
  B_2 = B_3*(w_c - w_0)^2 + (B_4/2)*(w_c - w_0);
  delta_w_p = (w_2 - w_1)/2;
  delta_w_s = w_s - w_c;
  w_sb1 = 2*delta_w_p - delta_w_s;
  %w_sb2 = 2*delta_w_p - delta_w_s;
  delta_w_sb = w_c - w_sb1;
  r = 2*gamma*sqrt(P_1*P_2);
dA_sb1dz = (1i/24)*A_sb1*(12*B_2*(delta_w_sb^2 - delta_w_p^2)... 
- 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_sb^4 - delta_w_p^4))+1i*gamma...
*P_1*(A_sb1 + conj(A_s))+ 1i*r*(conj(A_sb1) + A_s);
dA_sdz = (1i/24)*A_s*(12*B_2*(delta_w_s^2 - delta_w_p^2)... 
- 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_s^4 - delta_w_p^4))+ 1i*gamma...
*P_1*(conj(A_sb1) + A_s)+ 1i*r*(conj(A_1) + A_sb2);
dA_1dz = (1i/24)*A_1*(12*B_2*(delta_w_s^2 - delta_w_p^2)... 
+ 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_s^4 - delta_w_p^4))+1i*gamma...
*P_2*(conj(A_sb2) + A_1)+ 1i*r*(conj(A_s) + A_sb1);
dA_sb2dz = (1i/24)*A_sb2*(12*B_2*(delta_w_sb^2 - delta_w_p^2)... 
+ 4*B_3*(delta_w_sb^3-delta_w_p^3)+ B_4*(delta_w_sb^4 - delta_w_p^4))+1i*gamma...
*P_2*(conj(A_1) + A_sb2)+ 1i*r*(conj(A_sb1) + A_s);
     dAdz = [dA_sb1dz; dA_sdz; dA_1dz; dA_sb2dz];
  end
Related Question