Solving Second Order Differential Equation--Getting Nonsensical Answers!

Hi Everyone!

I am trying to solve the following equation using Matlab:

d^2(r)/dz^2 + (1/gamma)*dr/dz*dgamma/dz + (constant/gamma)*r = 0.

where gamma is a function of z but dgamma/dz is constant. I have set it up as shown below. However, my result appears to not depend on the dr/dz term. It gives the same result with or without that term! Any ideas as to what I am doing wrong?

Thanks so much!

~~~~~Code to establish function~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%===WRITE EQUATION AS FIRST ORDER SYSTEM===================================

%Let r1(z) be r(z)
%Let r2(z) be dr(z)/dz
%Let r3(z) be gamma
%Therefore r1prime(z)=dr(z)/dz=r2(z)

%Therefore our system of equations becomes
%==>r1prime = r' = r2
%==>r2prime = r'' = -(gamma^-1)*r2*dgamma-((w_p^2)/(2*gamma*c^2))*r1
%==>r3prime = r_3' = gamma-gamma_0

%===DEFINE AN M-FILE WITH THE SYSTEM OF EQUATIONS==========================

function rprime=M12_05_24_kBetaFunction(z,r)

%===USER INPUT=============================================================

plasma_density=3e19;
z_first_limit=500e-
lambda_0=800e-9; 

%Choose ramp profile
n_ramp=0; %Density step

%===DEFINE CONSTANTS=======================================================

e=1.6e-19; %Charge of electron [C].
c=3e8; %Speed of light [m/s].
epsilon=8.85e-12; %Permittivity of free space
m=9.11e-31; %Mass of electron [kg]

%===CALCULATE VALUES=======================================================

%Convert plasma density to SI units

plasma_density=plasma_density*100^3; %Convert plasma density to [m^-3].

%Set plasma density and energy gradient for each region

if z<z_first_limit
    n=plasma_density; %Set density inside plasma
    EnergyGradient=1.54e-20*sqrt(n); %Set energy gradient inside plasma
    %EnergyGradient=0;
else
    n=n_ramp; %Sets density on ramp
    EnergyGradient=0; %Assume no energy gain on ramp
end

%Define plasma frequencies
omega_0=2*pi*c/lambda_0;
omega_p=sqrt((e^2*n)/(m*epsilon));
omega_p0=sqrt((e^2*plasma_density)/(m*epsilon));

%Define gamma and energy at trapping
gamma_0=omega_0/omega_p0;
Energy_0=.511*(gamma_0-1)*1.6e-13; %[J]

%Set energy and gamma for each region

if z<z_first_limit
    Energy=Energy_0+EnergyGradient*z; %Calculate energy at each z [J]
    EnergyMeV=Energy/(1.6e-13); %Convert energy back to [MeV]
    gamma=EnergyMeV/.511+1; %Calcualte gamma at each z
else
    Energy=Energy_0+1.54e-20*sqrt(plasma_density)*z_first_limit; %Energy remains at the 
    %same value as when it leaves the plasma [J].
    EnergyMeV=Energy/(1.6e-13); %Convert energy back to [MeV]
    gamma=EnergyMeV/.511+1; %Calcualte gamma at each z
end

rprime=zeros(2,1);
rprime(1) = r(2);
rprime(2) = -(1/gamma)*r(2)*(1000/.511)*EnergyGradient-(omega_p^2)*r(1)/(2*c^2*gamma);

~~~~~Code to evaluate function~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%===CLOSE AND CLEAR EVERYTHING=============================================

close all
clear all
clc

%===USER INPUT=============================================================

lambda_0=800e-9; %Laser wavelength [m]
plasma_density=3e19; %Plasma density [cm^-3]

steps=100;
z_start=1e-6; %Location that electrons start oscillating [m]
z_exit=500e-6; %Location that plasma ends [m]
z_end=200e-6; %Length that electrons propagate in vacuum [m]

%Enter initial values of r1 and r2 (where again, r1=r(z) and r2=dr/dz).
r1_0=3e-6; %Distance off axis where electron bunch starts [m].
r2_0=0;

%===SET CONSTANTS==========================================================

e=1.6e-19; %Charge of electron [C]
c=3e8; %Speed of light [m/s].
epsilon=8.85e-12; %Permittivity of free space
m=9.11e-31; %Mass of electron [kg].

%===CALCULATE VALUES=======================================================

%Calculate plasma frequencies
omega_0=2*pi*c/lambda_0;
omega_p0=sqrt((e^2*plasma_density*100^3)/(m*epsilon));

%Define gamma at trapping
gamma_0=omega_0/omega_p0;

%Set initial values
r_0=[r1_0, r2_0]; %r_0(1) is r1_0 [m]. r_0(2) is r2_0 [dimenionsless]. 

ztotal=zeros(1,2*steps);
ztotal(1:steps)=linspace(z_start,z_exit,steps);
dz=(z_exit-z_end)/steps;
ztotal(steps+1:2*steps)=linspace(z_exit+dz,z_end+z_exit,steps);

%===ACTUALLY RUN ODE=======================================================

[Z,Result]=ode45(@M12_05_24_kBetaFunction,ztotal,r_0);

%Note: "Z" is all the values of z.  "Result(:,1)" is all the
%values of r.  "Results(:,2)" is all the values of dr/dz.

Answer 1

[Walter Roberson]

You indicate that dgamma/dz is constant . Is that the same constant as in your equation?

If they are then I find an analytic solution of

r(z) = C2 * BesselJ(0, 2*((constant*z + C1)/constant)^(1/2)) + 
         C3 * BesselY(0, 2*((constant*z + C1)/constant)^(1/2))
gamma(z) = constant*z + C1

In the above, C1, C2, C3 are constants of integration, and "constant" is the constant that dgamma(z)/dz is equal to.