clear all; close all

lam3 = 1e-6;
deff = 20e-12;
L=2e-2;
n=2.2;
c = 299792458;
e0 = 8.854187e-12;

phispan = [0 2*pi];

Ip = 100e7;
Is = 10e-2;

cmeqspan = [0 L];

ais = sqrt((2*Is)/(n*e0*c));
aip = sqrt((2*Ip)/(n*e0*c));

data = [];

for phi = linspace(0,2*pi)
    aphs = abs(ais)*exp(1i*phi);
    aInit = [aphs aip];
    [P,A] = ode45(@cmeq,cmeqspan,aInit);
    data = [data;A(0) A(end)];
end

function dA = cmeq(cmeqspan,A)
lam3 = 1e-6;
deff = 20e-12;
n=2.2;
c = 299792458;

w2 = 2*pi*c/lam3;
w1 = w2/2;
K1 = -(1i*w1*c*deff)/n;
K2 = -(1i*w2*c*deff)/n;

dA = zeros(2,1);
dA(1) = K1*A(2)*(A(1))';
dA(2) = K2*(A(1))^2;
end