system of differential equations

I have a system of differential equations as follows

where k1,k2... are constants. I tried solving it by using odeToVectorField and matlabFunction then ode45. However, I had "number of indeterminates exceeds the number of ODEs" error. Is there another way to solve this?

%% Photoionization Regime
syms P(z) R(z) ne(z) 

%Constants and initial conditions
l=8;
lambda0=775e-9;
n0=1;
n2=3e-23;
e=1.60217662e-19;                                     %elementary charge
Pcr=lambda0^2/(2*pi*n0*n2);                           %critical power
hbar=1.0545718e-34;                                   %h bar
epsilon0=8.85418782*1e-12;
mu0=4*pi*1e-7;
c=sqrt(1/epsilon0/mu0);                               %speed of light 
k0=2*pi/lambda0;                                      %wave number
w0=2*pi*c/lambda0;                                    %angular frequency
nn=2.7e25;                                            %neutral gas density
Uion=12.1*e;
% l=floor(Uion/hbar/w0+1);
sigmamp=6.4e-22;                                      %cross section
Imp=hbar*w0^2/sigmamp;                                %photoionization intensity
re=2.8e-15;                                            %classical electron radius
Pnl=lambda0^2/(2*pi*n2);                               %nonlinear focusing power
zf=2000;                                              %z final
intnum=5000;

%Initial conditions
alpha0=-.7;                                           %initial alpha
R0=1e-2;                                              %initial spot size
P0=Pcr;                                               %initial power
dRdz0=-2*alpha0/k0/R0; 
PR0=R0^2/P0;
I0=2*P0/pi/R0^2;                                      %initial intensity
dnedz0=2*pi*w0

%Differential equations
I=2*P/pi/R^2;
Wmp=2*pi*w0/factorial(l-1)*(I/Imp)^l;
dnedz=Wmp*nn/c;
dPdz=diff(P,1)==-c*pi/2*Uion*R^2/l*dnedz;
d2Rdz2=diff(R,2)==4/k0^2/R^3*(1-P/Pnl+2*pi*l/(l+1)^2*...
    re*R^2*ne)-(l-1)/2/l/R^3*(R^2*diff(R^2/P*dPdz,1)+(l-1)/2/l...
    *(R^2/P*dPdz)^2);
eqn=[d2Rdz2 dPdz dnedz];
v=odeToVectorField(eqn);     %HERE I GET THE ERROR

[darova]

Please attach the code

[Ayberk Cakmak]

I resubmitted with the code attached. Thanks in advance.

Answer 1

[darova]

Here is the problem: there is second derivative of , but there is no expression for

Looks like 4 variables, but only 3 equations