I am developing a code which uses NDSolve to find solutions to two sets of two simultaneous differential equations, I would like to use a Do loop or a For loop to change some variables and repeat the calculations, to find the point where the sets diverge, however if I use a Do loop or a For loop, even without changing any other information, for instance, using a Do loop to do the calculations only once, NDSolve no longer works.

m = 9.1*10^-28; (*Mass of Electron in g*)(*Initialising constants, \

starting with universal constants*)

c = 3*10^10;(*Speed of light in cm/s*)

e = 4.8*10^-10;(*Electron charge in StatColoumb*)

\[Lambda] = 10^-4;(*Wavelength of laser in cm*)

r0 = 30*10^-4;(*Spot size in cm*)

n0 = 6*10^17;(*Electron density in cm^-3*)

energy = 10^5;(*Energy of laser pulse in erg*)

FWHM = 100*10^-15;(*Full-width half maximum of power distribution in \

seconds*)

\[Tau] = FWHM/(2 Sqrt[Log[2]]) // N;(*Time Constant in seconds*)

\[Omega]0 = 2*Pi*c/\[Lambda];

\[Omega]p = Sqrt[4*Pi*(e^2)*n0/m];(*Plasma Frequency in rad.s^-1*)

kp = \[Omega]p/c // N;

\[Tau]p = 2*Pi/\[Omega]p;

vp = c (1 - (\[Omega]p/\[Omega]0)^2)^0.5;

kR = r0*kp;

kZ = kp*vp*\[Tau];

P0 = (energy/FWHM)*Sqrt[2*Log[2]/Pi];(*Peak Power in erg/s*)

I0 = 2 P0/(Pi*r0^2);(*Peak Intensity in erg/(cm^2/s)*)

a0 = Sqrt[(2/(c*Pi))*(\[Lambda]*e/(m*c^2))^2*I0];(*Amplitude constant*)

asquared1[Z_, R_, T_] :=

a0^2*Exp[-(((Z - T)/kZ)^2)]*Exp[-2 ((R/kR)^2)];

Tmin = -20;

Tmax = 20;

Zmin = 0;

Zmax = 0.1;

Rmin = -0.1;

Rmax = 0.1;

T0 = -20;

Do[

s = NDSolve[{D[

u1[Z, R, T], {T, 2}]*(1 + 1.5*D[u1[Z, R, T], T]^2 +

0.5*D[v1[Z, R, T], T]^2) + u1[Z, R, T] -

D[v1[Z, R, T], T]*D[u1[Z, R, T], T]*v1[Z, R, T] == -0.5*

D[asquared1[Z, R, T], R],

D[v1[Z, R, T], {T, 2}]*(1 + 1.5*D[v1[Z, R, T], T]^2 +

0.5*D[u1[Z, R, T], T]^2) + v1[Z, R, T] -

D[u1[Z, R, T], T]*D[v1[Z, R, T], T]*u1[Z, R, T] == -0.5*

D[asquared1[Z, R, T], Z],

(D[u1[Z, R, T], T] /. T -> T0) == 0, v1[Z, R, T0] == 0,

u1[Z, R, T0] == 0, (D[v1[Z, R, T], T] /. T -> T0) == 0}, {u1,

v1}, {Z, Zmin, Zmax}, {R, Rmin, Rmax}, {T, Tmin, Tmax},

MaxSteps -> Infinity, MaxStepSize -> {0.1, 0.1, 0.1},

Method -> "ExplicitRungeKutta"]

First[NMaximize[{v1[0, 0, T] /. s[[1]], 10 < T < 20}, T]]

, {1}]

Removing the Do loop, the code works, but if I add it or a For loop, it gives the errros NMaximize::nnum: and Thread::tdlen: Could anyone suggest where I'm going wrong? Thank you.

Sorry for any mistakes with formatting, this is my first time using this website.