 CoCalc Public FilesMIWsim1DjacobianV1b.sagews
Authors: David Cyganski, Bill Page
Views : 56
Description: Implementation of Wiseman Many-Interacting-Worlds simulation for 1D, 1 Particle systems
Compute Environment: Ubuntu 18.04 (Deprecated)

###MIWsim1D###

• using lsode with Jacobian

Author: D Cyganski and Bill Page

Created: 2015-05-17

Revised: 2015-06-08

####Implementation of Wiseman Many-Interacting-Worlds simulation for 1D, 1 Particle systems####

• single particle with gaussian distribution and nonzero mean velocity, colliding with a potential barrier

The future track of the many worlds position state $Xp(t)=[x_1(t)...x_{N*J*D}(t)]$ is determined in general in MIW by system of $N*J*D$ coupled non-linear second order partial differential equations, where $N$=number of worlds, $J$=number of particles, $D$=number of dimensions. $m_k*\ddot{x}_{kn}(t) = f_k(x_n(t)) + r_{kN}(x_n(t);Xp(t))$ where $k$ indexes the particle and $n$ the world.

Such a system of $N*J*D$ second order non-linear PDEs can always be written as a system of of $2*N*J*D$ first order linear PDEs:

$\dot{x}_{kn}(t) = y_{kn}(t) \\ \dot{y}_{kn}(t) = (1/m_k)*(f_k(x_n(t)) + r_{kN}(x_n(t))$

For our 1 D case, there will be $2N$ equations.

Setup Sage environment

%sage
# Better version of FriCAS/Axiom mode for SageMathCloud
os.environ['PATH'] = '%s/bin:%s'%(os.environ['HOME'],os.environ['PATH'])
from fricas import fricas
execfile('fricas_md.py')
%default_mode fricas_md

)version
)clear completely
)set output tex on
)set output algebra off
)set message type off


Value = "FriCAS 2014-12-18 compiled at Tue May 5 18:13:01 UTC 2015"

All user variables and function definitions have been cleared.

All )browse facility databases have been cleared.

Internally cached functions and constructors have been cleared.

)clear completely is finished.

Classical Force centered at $x=c$ that repells the oncoming particles

--F(x) == b*(x-c)/(a^3*%pi)*exp(-(x-c)^2/a^2);
F(x) == 9*sech(5*(x-c))^2*tanh(5*(x-c))
F(x)


Compiling function F with type Variable(x) -> Expression(Integer)

$$9 \ {{{\mathrm{sech} \left( {{{5 \ x} -{5 \ c}}} \right)}} ^{2}} \ {\mathrm{tanh} \left( {{{5 \ x} -{5 \ c}}} \right)}$$

complexNormalize htrigs F(x)


$${{{36} \ {{{{e} ^{{{5 \ x} -{5 \ c}}}}} ^{4}}} -{{36} \ {{{{e} ^{{{5 \ x} -{5 \ c}}}}} ^{2}}}} \over {{{{{e} ^{{{5 \ x} -{5 \ c}}}}} ^{6}}+{3 \ {{{{e} ^{{{5 \ x} -{5 \ c}}}}} ^{4}}}+{3 \ {{{{e} ^{{{5 \ x} -{5 \ c}}}}} ^{2}}}+1}$$

%sage
plot(fricas("eval(F(x),[a=0.25,c=5,b=0.7])").sage(),(x,0,10),color="red") Potential

%sage
plot(fricas("eval(-integrate(F(x),x),[a=0.25,c=5,b=0.7])").sage(),(x,0,10),color="red") Position of particle in the $i^{th}$ world: $x(i)$

x:=operator('x);


Many Worlds Interaction

Equations 24,25 of HEW paper

R(x,i) == (hbar^2/(4*m))*(sigma(x,i+1)-sigma(x,i))
sigma(x,i) == sif(x,i)^2*(sif(x,i+1)-2*sif(x,i)+sif(x,i-1))
sif(x,i) == 1/(x(i)-x(i-1))


Equations of Motion

xdd(i)==(F(x(i))+R(x,i))/m


The non-constant elements of the Jacobian matrix are formed from derivates with respect to x

jac1(i,j)==realElementary eval(D(eval(xdd(i),x(i+j)=%x),%x),%x=x(i+j))


Sage routine to convert FriCAS expression to Octave function

%sage
# fold lines for octave
def split_lines(fr):
i=0; cpl=80; f=""
while i<len(fr):
if i+cpl>=len(fr):
f+= "    "+fr[i:i+cpl]+";\n"
i=i+cpl
else:
line=fr[i:i+cpl]
p=line.rfind('*', 0, cpl)+1
if p==0:
p=line.rindex('+', 0, cpl)+1
f+= "    "+fr[i:i+p]+"...\n"
i=i+p
return f[4:len(f)-1]

%sage
# name of function defined by fricas expression ex
def fricas2octave(name,ex):
fr=fricas(ex).unparsed_input_form()
# cleanup fricas expresssion
fr = fr.replace("\n","")
fr = fr.replace("\r","")
fr = fr.replace(" -","-")
fr = fr.replace("+-","-")
# boundry conditions
fr = fr.replace("x(i-2)","xl2")
fr = fr.replace("x(i-1)","xl1")
fr = fr.replace("x(i+1)","xr1")
fr = fr.replace("x(i+2)","xr2")
f = """\
function result=%s(x,i)
global l2 l1 p1 p2;
global N m hbar a b c;
if (i-2<1)
xl2=l2;
else
xl2=x(i-2);
endif;
if (i-1<1)
xl1=l1;
else
xl1=x(i-1);
endif;
if (i+2>N)
xr2=p2;
else
xr2=x(i+2);
endif;
if (i+1>N)
xr1=p1;
else
xr1=x(i+1);
endif;
result=%s
endfunction
"""%(name,split_lines(fr))
ff = open('%s.m'%name, 'w')
ff.write(f)
ff.close()
return f

%sage
# Use a better version of the octave interface
from octave import octave
%default_mode octave


Generate initial distribution of particle states, including position,x (X0), and velocity, xdot (Y0)

clear all;
version

ans = 3.8.2

For a simple test, distribute 101 values equally over a window and set all velocities to zero

#xvalues=[-1:.05:1]'; #Prob2
#xvalues=[-1:.025:1]'; #Prob3
#xvalues=[-0.5:.025:0.5]'; #Prob4


Use inverse error function to generate a Gaussian distribution of trajectories

#xuniform=[-0.95:0.1:0.95]'; %For the simple test
#xuniform=[-0.95:0.022|:0.95]';
xuniform=[-0.99:0.15:0.99]';
#xvalues=erfcinv(xuniform);
xvalues=erfinv(xuniform);
global N; %Number of worlds
N = size(xvalues,1)
X0 = xvalues;
#Y0 = zeros(N,1); %zero initial velocities
#Y0 = normrnd(1, 0.2, N, 1);
Y0 = ones(N,1); %non-zero initial velocities

N = 14

Plot the input trajectory density

Pin=(N*(xvalues(2:N,1)-xvalues(1:(N-1),1))).^(-1);

%sage list_plot(matrix(RR,octave("[xvalues(2:N),Pin]")),plotjoined=true) Constants

global m hbar;
m = 1; %mass of particle
hbar = 0.1; %Planck constant %.02 works but nonguassian after 500 iterations
global a b c;
a = 0.25;
b = 0.7;
c = 5;


Boundry conditions

global l2 l1 p1 p2;
l2 = -100000;
l1 =  -10000;
p1 =   10000;
p2 =  100000;

%sage
octave.eval("global F = @(x) "+fricas("F('x)").unparsed_input_form().replace("pi()","pi"))

''

Plot the force

Xrange = 0:0.05:10;

%sage list_plot(matrix(RR,octave("[Xrange',arrayfun(F,Xrange')]")),plotjoined=true) %sage print fricas2octave('R2','R(x,i)')

function result=R2(x,i) global l2 l1 p1 p2; global N m hbar a b c; if (i-2<1) xl2=l2; else xl2=x(i-2); endif; if (i-1<1) xl1=l1; else xl1=x(i-1); endif; if (i+2>N) xr2=p2; else xr2=x(i+2); endif; if (i+1>N) xr1=p1; else xr1=x(i+1); endif; result=(((-1*hbar^2*x(i)+(3*hbar^2*xl1-2*hbar^2*xl2))*xr1^3+(3*hbar^2*x(i)^2+(-10*... hbar^2*xl1+7*hbar^2*xl2)*x(i)+(hbar^2*xl1^2-1*hbar^2*xl2*xl1))*xr1^2+(-3*hbar^2*... x(i)^3+(12*hbar^2*xl1-9*hbar^2*xl2)*x(i)^2+(-4*hbar^2*xl1^2+4*hbar^2*xl2*xl1)*... x(i)+(hbar^2*xl1^3-1*hbar^2*xl2*xl1^2))*xr1+(hbar^2*x(i)^4+(-7*hbar^2*xl1+6*... hbar^2*xl2)*x(i)^3+(9*hbar^2*xl1^2-9*hbar^2*xl2*xl1)*x(i)^2+(-7*hbar^2*xl1^3+7*... hbar^2*xl2*xl1^2)*x(i)+(2*hbar^2*xl1^4-2*hbar^2*xl2*xl1^3)))*xr2+((hbar^2*... x(i)+(-3*hbar^2*xl1+2*hbar^2*xl2))*xr1^4+(-3*hbar^2*x(i)^2+(10*hbar^2*xl1-7*... hbar^2*xl2)*x(i)+(-1*hbar^2*xl1^2+hbar^2*xl2*xl1))*xr1^3+(3*hbar^2*x(i)^3+(-12*... hbar^2*xl1+9*hbar^2*xl2)*x(i)^2+(4*hbar^2*xl1^2-4*hbar^2*xl2*xl1)*x(i)+(-1*... hbar^2*xl1^3+hbar^2*xl2*xl1^2))*xr1^2+(-1*hbar^2*x(i)^4+(8*hbar^2*xl1-7*hbar^2*... xl2)*x(i)^3+(-12*hbar^2*xl1^2+12*hbar^2*xl2*xl1)*x(i)^2+(10*hbar^2*xl1^3-10*... hbar^2*xl2*xl1^2)*x(i)+(-3*hbar^2*xl1^4+3*hbar^2*xl2*xl1^3))*xr1+((-1*hbar^2*... xl1+hbar^2*xl2)*x(i)^4+(3*hbar^2*xl1^2-3*hbar^2*xl2*xl1)*x(i)^3+(-3*hbar^2*... xl1^3+3*hbar^2*xl2*xl1^2)*x(i)^2+(hbar^2*xl1^4-1*hbar^2*xl2*xl1^3)*... x(i))))/((((4*m*xl1-4*m*xl2)*x(i)^3+(-12*m*xl1^2+12*m*xl2*xl1)*x(i)^2+(12*m*... xl1^3-12*m*xl2*xl1^2)*x(i)+(-4*m*xl1^4+4*m*xl2*xl1^3))*xr1^3+((-12*m*xl1+12*m*... xl2)*x(i)^4+(36*m*xl1^2-36*m*xl2*xl1)*x(i)^3+(-36*m*xl1^3+36*m*xl2*xl1^2)*... x(i)^2+(12*m*xl1^4-12*m*xl2*xl1^3)*x(i))*xr1^2+((12*m*xl1-12*m*xl2)*x(i)^5+(-36*... m*xl1^2+36*m*xl2*xl1)*x(i)^4+(36*m*xl1^3-36*m*xl2*xl1^2)*x(i)^3+(-12*m*xl1^4+12*... m*xl2*xl1^3)*x(i)^2)*xr1+((-4*m*xl1+4*m*xl2)*x(i)^6+(12*m*xl1^2-12*m*xl2*xl1)*... x(i)^5+(-12*m*xl1^3+12*m*xl2*xl1^2)*x(i)^4+(4*m*xl1^4-4*m*xl2*xl1^3)*x(i)^3))*... xr2+(((-4*m*xl1+4*m*xl2)*x(i)^3+(12*m*xl1^2-12*m*xl2*xl1)*x(i)^2+(-12*m*... xl1^3+12*m*xl2*xl1^2)*x(i)+(4*m*xl1^4-4*m*xl2*xl1^3))*xr1^4+((12*m*xl1-12*m*... xl2)*x(i)^4+(-36*m*xl1^2+36*m*xl2*xl1)*x(i)^3+(36*m*xl1^3-36*m*xl2*xl1^2)*... x(i)^2+(-12*m*xl1^4+12*m*xl2*xl1^3)*x(i))*xr1^3+((-12*m*xl1+12*m*xl2)*... x(i)^5+(36*m*xl1^2-36*m*xl2*xl1)*x(i)^4+(-36*m*xl1^3+36*m*xl2*xl1^2)*x(i)^3+(12*... m*xl1^4-12*m*xl2*xl1^3)*x(i)^2)*xr1^2+((4*m*xl1-4*m*xl2)*x(i)^6+(-12*m*xl1^2+12*... m*xl2*xl1)*x(i)^5+(12*m*xl1^3-12*m*xl2*xl1^2)*x(i)^4+(-4*m*xl1^4+4*m*xl2*xl1^3)*... x(i)^3)*xr1)); endfunction
R2(X0,3)

ans = 0.00280219
function r = R(x,i)
global m hbar;
r = (hbar^2/(4*m))*(sigma(x,i+1)-sigma(x,i));
endfunction;

function s=sigma(x,i)
s = sif(x,i)^2*(sif(x,i+1)-2*sif(x,i)+sif(x,i-1));
endfunction;

function den = sif(x,i)
global N
if (i>N)||(i<2)
den=0;
else
den=1/(x(i)-x(i-1));
endif
endfunction;

R2(X0,3)

ans = 0.00280219
%sage
X=map(eval,str(octave('X0')).splitlines()[2:-1]);X

[-1.82139, -0.993536, -0.71787, -0.522444, -0.360676, -0.216008, -0.0799303, 0.0532238, 0.18831, 0.330713, 0.488122, 0.674697, 0.926719, 1.45222]
%sage
print fricas.eval("X:=%s;"%(fricas(X).name()))

 
%fricas_md
eval(R(X,3),[hbar=0.1,m=1])


Compiling function sif with type (List(Float),PositiveInteger) -> Float

Compiling function sif with type (List(Float),Integer) -> Float

Compiling function sigma with type (List(Float),PositiveInteger) -> Float

Compiling function R with type (List(Float),PositiveInteger) -> Fraction(Polynomial(Float))

$$0.0027954071_1864352625_97$$

Equations of motion in numeric form

function xdot = deq(x,t)
global F N;
for i=1:N
xdot(i) = x(N+i);
xdot(N+i) = F(x(i))+R(x,i);
endfor
endfunction


Functions to compute the Jacobian matrix

%sage print fricas2octave('Jm2','jac1(i,-2)')

function result=Jm2(x,i) global l2 l1 p1 p2; global N m hbar a b c; if (i-2<1) xl2=l2; else xl2=x(i-2); endif; if (i-1<1) xl1=l1; else xl1=x(i-1); endif; if (i+2>N) xr2=p2; else xr2=x(i+2); endif; if (i+1>N) xr1=p1; else xr1=x(i+1); endif; result=(-1*hbar^2)/((4*m^2*xl1^2-8*m^2*xl2*xl1+4*m^2*xl2^2)*x(i)^2+(-8*m^2*xl1^3+16*... m^2*xl2*xl1^2-8*m^2*xl2^2*xl1)*x(i)+(4*m^2*xl1^4-8*m^2*xl2*xl1^3+4*m^2*xl2^2*... xl1^2)); endfunction
%sage print fricas2octave('Jm1','jac1(i,-1)')

function result=Jm1(x,i) global l2 l1 p1 p2; global N m hbar a b c; if (i-2<1) xl2=l2; else xl2=x(i-2); endif; if (i-1<1) xl1=l1; else xl1=x(i-1); endif; if (i+2>N) xr2=p2; else xr2=x(i+2); endif; if (i+1>N) xr1=p1; else xr1=x(i+1); endif; result=((hbar^2*x(i)^2+(-4*hbar^2*xl1+2*hbar^2*xl2)*x(i)+(9*hbar^2*xl1^2-14*hbar^2*xl2*... xl1+6*hbar^2*xl2^2))*xr1^2+(-2*hbar^2*x(i)^3+(8*hbar^2*xl1-4*hbar^2*xl2)*... x(i)^2+(-20*hbar^2*xl1^2+32*hbar^2*xl2*xl1-14*hbar^2*xl2^2)*x(i)+(2*hbar^2*... xl1^3-4*hbar^2*xl2*xl1^2+2*hbar^2*xl2^2*xl1))*xr1+(hbar^2*x(i)^4+(-4*hbar^2*... xl1+2*hbar^2*xl2)*x(i)^3+(12*hbar^2*xl1^2-20*hbar^2*xl2*xl1+9*hbar^2*xl2^2)*... x(i)^2+(-4*hbar^2*xl1^3+8*hbar^2*xl2*xl1^2-4*hbar^2*xl2^2*xl1)*x(i)+(hbar^2*... xl1^4-2*hbar^2*xl2*xl1^3+hbar^2*xl2^2*xl1^2)))/(((4*m^2*xl1^2-8*m^2*xl2*xl1+4*... m^2*xl2^2)*x(i)^4+(-16*m^2*xl1^3+32*m^2*xl2*xl1^2-16*m^2*xl2^2*xl1)*x(i)^3+(24*... m^2*xl1^4-48*m^2*xl2*xl1^3+24*m^2*xl2^2*xl1^2)*x(i)^2+(-16*m^2*xl1^5+32*m^2*xl2*... xl1^4-16*m^2*xl2^2*xl1^3)*x(i)+(4*m^2*xl1^6-8*m^2*xl2*xl1^5+4*m^2*xl2^2*xl1^4))*... xr1^2+((-8*m^2*xl1^2+16*m^2*xl2*xl1-8*m^2*xl2^2)*x(i)^5+(32*m^2*xl1^3-64*m^2*... xl2*xl1^2+32*m^2*xl2^2*xl1)*x(i)^4+(-48*m^2*xl1^4+96*m^2*xl2*xl1^3-48*m^2*xl2^2*... xl1^2)*x(i)^3+(32*m^2*xl1^5-64*m^2*xl2*xl1^4+32*m^2*xl2^2*xl1^3)*x(i)^2+(-8*m^2*... xl1^6+16*m^2*xl2*xl1^5-8*m^2*xl2^2*xl1^4)*x(i))*xr1+((4*m^2*xl1^2-8*m^2*xl2*... xl1+4*m^2*xl2^2)*x(i)^6+(-16*m^2*xl1^3+32*m^2*xl2*xl1^2-16*m^2*xl2^2*xl1)*... x(i)^5+(24*m^2*xl1^4-48*m^2*xl2*xl1^3+24*m^2*xl2^2*xl1^2)*x(i)^4+(-16*m^2*... xl1^5+32*m^2*xl2*xl1^4-16*m^2*xl2^2*xl1^3)*x(i)^3+(4*m^2*xl1^6-8*m^2*xl2*... xl1^5+4*m^2*xl2^2*xl1^4)*x(i)^2)); endfunction
%sage print fricas2octave('J0','jac1(i,0)')

function result=J0(x,i) global l2 l1 p1 p2; global N m hbar a b c; if (i-2<1) xl2=l2; else xl2=x(i-2); endif; if (i-1<1) xl1=l1; else xl1=x(i-1); endif; if (i+2>N) xr2=p2; else xr2=x(i+2); endif; if (i+1>N) xr1=p1; else xr1=x(i+1); endif; result=((((hbar^2*x(i)+(-4*hbar^2*xl1+3*hbar^2*xl2))*xr1^4+(-4*hbar^2*x(i)^2+(17*... hbar^2*xl1-13*hbar^2*xl2)*x(i)+(-1*hbar^2*xl1^2+hbar^2*xl2*xl1))*xr1^3+(6*... hbar^2*x(i)^3+(-28*hbar^2*xl1+22*hbar^2*xl2)*x(i)^2+(5*hbar^2*xl1^2-5*hbar^2*... xl2*xl1)*x(i)+(-1*hbar^2*xl1^3+hbar^2*xl2*xl1^2))*xr1^2+(-4*hbar^2*x(i)^4+(22*... hbar^2*xl1-18*hbar^2*xl2)*x(i)^3+(-10*hbar^2*xl1^2+10*hbar^2*xl2*xl1)*x(i)^2+(5*... hbar^2*xl1^3-5*hbar^2*xl2*xl1^2)*x(i)+(-1*hbar^2*xl1^4+hbar^2*xl2*xl1^3))*... xr1+(hbar^2*x(i)^5+(-10*hbar^2*xl1+9*hbar^2*xl2)*x(i)^4+(18*hbar^2*xl1^2-18*... hbar^2*xl2*xl1)*x(i)^3+(-22*hbar^2*xl1^3+22*hbar^2*xl2*xl1^2)*x(i)^2+(13*hbar^2*... xl1^4-13*hbar^2*xl2*xl1^3)*x(i)+(-3*hbar^2*xl1^5+3*hbar^2*xl2*xl1^4)))*xr2+((-1*... hbar^2*x(i)+(4*hbar^2*xl1-3*hbar^2*xl2))*xr1^5+(4*hbar^2*x(i)^2+(-17*hbar^2*... xl1+13*hbar^2*xl2)*x(i)+(hbar^2*xl1^2-1*hbar^2*xl2*xl1))*xr1^4+(-6*hbar^2*... x(i)^3+(28*hbar^2*xl1-22*hbar^2*xl2)*x(i)^2+(-5*hbar^2*xl1^2+5*hbar^2*xl2*xl1)*... x(i)+(hbar^2*xl1^3-1*hbar^2*xl2*xl1^2))*xr1^3+(4*hbar^2*x(i)^4+(-22*hbar^2*... xl1+18*hbar^2*xl2)*x(i)^3+(10*hbar^2*xl1^2-10*hbar^2*xl2*xl1)*x(i)^2+(-5*hbar^2*... xl1^3+5*hbar^2*xl2*xl1^2)*x(i)+(hbar^2*xl1^4-1*hbar^2*xl2*xl1^3))*xr1^2+(-1*... hbar^2*x(i)^5+(11*hbar^2*xl1-10*hbar^2*xl2)*x(i)^4+(-22*hbar^2*xl1^2+22*hbar^2*... xl2*xl1)*x(i)^3+(28*hbar^2*xl1^3-28*hbar^2*xl2*xl1^2)*x(i)^2+(-17*hbar^2*... xl1^4+17*hbar^2*xl2*xl1^3)*x(i)+(4*hbar^2*xl1^5-4*hbar^2*xl2*xl1^4))*xr1+((-1*... hbar^2*xl1+hbar^2*xl2)*x(i)^5+(4*hbar^2*xl1^2-4*hbar^2*xl2*xl1)*x(i)^4+(-6*... hbar^2*xl1^3+6*hbar^2*xl2*xl1^2)*x(i)^3+(4*hbar^2*xl1^4-4*hbar^2*xl2*xl1^3)*... x(i)^2+(-1*hbar^2*xl1^5+hbar^2*xl2*xl1^4)*x(i))))*exp(5*x(i)-5*c)^8+((((-720*m*... xl1+720*m*xl2)*x(i)^4+(2880*m*xl1^2-2880*m*xl2*xl1)*x(i)^3+(-4320*m*xl1^3+4320*... m*xl2*xl1^2)*x(i)^2+(2880*m*xl1^4-2880*m*xl2*xl1^3+4*hbar^2)*x(i)+(-720*m*... xl1^5+720*m*xl2*xl1^4-16*hbar^2*xl1+12*hbar^2*xl2))*xr1^4+((2880*m*xl1-2880*m*... xl2)*x(i)^5+(-11520*m*xl1^2+11520*m*xl2*xl1)*x(i)^4+(17280*m*xl1^3-17280*m*xl2*... xl1^2)*x(i)^3+(-11520*m*xl1^4+11520*m*xl2*xl1^3-16*hbar^2)*x(i)^2+(2880*m*... xl1^5-2880*m*xl2*xl1^4+68*hbar^2*xl1-52*hbar^2*xl2)*x(i)+(-4*hbar^2*xl1^2+4*... hbar^2*xl2*xl1))*xr1^3+((-4320*m*xl1+4320*m*xl2)*x(i)^6+(17280*m*xl1^2-17280*m*... xl2*xl1)*x(i)^5+(-25920*m*xl1^3+25920*m*xl2*xl1^2)*x(i)^4+(17280*m*xl1^4-17280*... m*xl2*xl1^3+24*hbar^2)*x(i)^3+(-4320*m*xl1^5+4320*m*xl2*xl1^4-112*hbar^2*xl1+88*... hbar^2*xl2)*x(i)^2+(20*hbar^2*xl1^2-20*hbar^2*xl2*xl1)*x(i)+(-4*hbar^2*xl1^3+4*... hbar^2*xl2*xl1^2))*xr1^2+((2880*m*xl1-2880*m*xl2)*x(i)^7+(-11520*m*xl1^2+11520*... m*xl2*xl1)*x(i)^6+(17280*m*xl1^3-17280*m*xl2*xl1^2)*x(i)^5+(-11520*m*... xl1^4+11520*m*xl2*xl1^3-16*hbar^2)*x(i)^4+(2880*m*xl1^5-2880*m*xl2*xl1^4+88*... hbar^2*xl1-72*hbar^2*xl2)*x(i)^3+(-40*hbar^2*xl1^2+40*hbar^2*xl2*xl1)*... x(i)^2+(20*hbar^2*xl1^3-20*hbar^2*xl2*xl1^2)*x(i)+(-4*hbar^2*xl1^4+4*hbar^2*xl2*... xl1^3))*xr1+((-720*m*xl1+720*m*xl2)*x(i)^8+(2880*m*xl1^2-2880*m*xl2*xl1)*... x(i)^7+(-4320*m*xl1^3+4320*m*xl2*xl1^2)*x(i)^6+(2880*m*xl1^4-2880*m*xl2*xl1^3+4*... hbar^2)*x(i)^5+(-720*m*xl1^5+720*m*xl2*xl1^4-40*hbar^2*xl1+36*hbar^2*xl2)*... x(i)^4+(72*hbar^2*xl1^2-72*hbar^2*xl2*xl1)*x(i)^3+(-88*hbar^2*xl1^3+88*hbar^2*... xl2*xl1^2)*x(i)^2+(52*hbar^2*xl1^4-52*hbar^2*xl2*xl1^3)*x(i)+(-12*hbar^2*... xl1^5+12*hbar^2*xl2*xl1^4)))*xr2+(((720*m*xl1-720*m*xl2)*x(i)^4+(-2880*m*... xl1^2+2880*m*xl2*xl1)*x(i)^3+(4320*m*xl1^3-4320*m*xl2*xl1^2)*x(i)^2+(-2880*m*... xl1^4+2880*m*xl2*xl1^3-4*hbar^2)*x(i)+(720*m*xl1^5-720*m*xl2*xl1^4+16*hbar^2*... xl1-12*hbar^2*xl2))*xr1^5+((-2880*m*xl1+2880*m*xl2)*x(i)^5+(11520*m*xl1^2-11520*... m*xl2*xl1)*x(i)^4+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^3+(11520*m*... xl1^4-11520*m*xl2*xl1^3+16*hbar^2)*x(i)^2+(-2880*m*xl1^5+2880*m*xl2*xl1^4-68*... hbar^2*xl1+52*hbar^2*xl2)*x(i)+(4*hbar^2*xl1^2-4*hbar^2*xl2*xl1))*xr1^4+((4320*... m*xl1-4320*m*xl2)*x(i)^6+(-17280*m*xl1^2+17280*m*xl2*xl1)*x(i)^5+(25920*m*... xl1^3-25920*m*xl2*xl1^2)*x(i)^4+(-17280*m*xl1^4+17280*m*xl2*xl1^3-24*hbar^2)*... x(i)^3+(4320*m*xl1^5-4320*m*xl2*xl1^4+112*hbar^2*xl1-88*hbar^2*xl2)*x(i)^2+(-20*... hbar^2*xl1^2+20*hbar^2*xl2*xl1)*x(i)+(4*hbar^2*xl1^3-4*hbar^2*xl2*xl1^2))*... xr1^3+((-2880*m*xl1+2880*m*xl2)*x(i)^7+(11520*m*xl1^2-11520*m*xl2*xl1)*... x(i)^6+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^5+(11520*m*xl1^4-11520*m*xl2*... xl1^3+16*hbar^2)*x(i)^4+(-2880*m*xl1^5+2880*m*xl2*xl1^4-88*hbar^2*xl1+72*hbar^2*... xl2)*x(i)^3+(40*hbar^2*xl1^2-40*hbar^2*xl2*xl1)*x(i)^2+(-20*hbar^2*xl1^3+20*... hbar^2*xl2*xl1^2)*x(i)+(4*hbar^2*xl1^4-4*hbar^2*xl2*xl1^3))*xr1^2+((720*m*... xl1-720*m*xl2)*x(i)^8+(-2880*m*xl1^2+2880*m*xl2*xl1)*x(i)^7+(4320*m*xl1^3-4320*... m*xl2*xl1^2)*x(i)^6+(-2880*m*xl1^4+2880*m*xl2*xl1^3-4*hbar^2)*x(i)^5+(720*m*... xl1^5-720*m*xl2*xl1^4+44*hbar^2*xl1-40*hbar^2*xl2)*x(i)^4+(-88*hbar^2*xl1^2+88*... hbar^2*xl2*xl1)*x(i)^3+(112*hbar^2*xl1^3-112*hbar^2*xl2*xl1^2)*x(i)^2+(-68*... hbar^2*xl1^4+68*hbar^2*xl2*xl1^3)*x(i)+(16*hbar^2*xl1^5-16*hbar^2*xl2*xl1^4))*... xr1+((-4*hbar^2*xl1+4*hbar^2*xl2)*x(i)^5+(16*hbar^2*xl1^2-16*hbar^2*xl2*xl1)*... x(i)^4+(-24*hbar^2*xl1^3+24*hbar^2*xl2*xl1^2)*x(i)^3+(16*hbar^2*xl1^4-16*hbar^2*... xl2*xl1^3)*x(i)^2+(-4*hbar^2*xl1^5+4*hbar^2*xl2*xl1^4)*x(i))))*exp(5*x(i)-5*... c)^6+((((2880*m*xl1-2880*m*xl2)*x(i)^4+(-11520*m*xl1^2+11520*m*xl2*xl1)*... x(i)^3+(17280*m*xl1^3-17280*m*xl2*xl1^2)*x(i)^2+(-11520*m*xl1^4+11520*m*xl2*... xl1^3+6*hbar^2)*x(i)+(2880*m*xl1^5-2880*m*xl2*xl1^4-24*hbar^2*xl1+18*hbar^2*... xl2))*xr1^4+((-11520*m*xl1+11520*m*xl2)*x(i)^5+(46080*m*xl1^2-46080*m*xl2*xl1)*... x(i)^4+(-69120*m*xl1^3+69120*m*xl2*xl1^2)*x(i)^3+(46080*m*xl1^4-46080*m*xl2*... xl1^3-24*hbar^2)*x(i)^2+(-11520*m*xl1^5+11520*m*xl2*xl1^4+102*hbar^2*xl1-78*... hbar^2*xl2)*x(i)+(-6*hbar^2*xl1^2+6*hbar^2*xl2*xl1))*xr1^3+((17280*m*xl1-17280*... m*xl2)*x(i)^6+(-69120*m*xl1^2+69120*m*xl2*xl1)*x(i)^5+(103680*m*xl1^3-103680*m*... xl2*xl1^2)*x(i)^4+(-69120*m*xl1^4+69120*m*xl2*xl1^3+36*hbar^2)*x(i)^3+(17280*m*... xl1^5-17280*m*xl2*xl1^4-168*hbar^2*xl1+132*hbar^2*xl2)*x(i)^2+(30*hbar^2*... xl1^2-30*hbar^2*xl2*xl1)*x(i)+(-6*hbar^2*xl1^3+6*hbar^2*xl2*xl1^2))*... xr1^2+((-11520*m*xl1+11520*m*xl2)*x(i)^7+(46080*m*xl1^2-46080*m*xl2*xl1)*... x(i)^6+(-69120*m*xl1^3+69120*m*xl2*xl1^2)*x(i)^5+(46080*m*xl1^4-46080*m*xl2*... xl1^3-24*hbar^2)*x(i)^4+(-11520*m*xl1^5+11520*m*xl2*xl1^4+132*hbar^2*xl1-108*... hbar^2*xl2)*x(i)^3+(-60*hbar^2*xl1^2+60*hbar^2*xl2*xl1)*x(i)^2+(30*hbar^2*... xl1^3-30*hbar^2*xl2*xl1^2)*x(i)+(-6*hbar^2*xl1^4+6*hbar^2*xl2*xl1^3))*... xr1+((2880*m*xl1-2880*m*xl2)*x(i)^8+(-11520*m*xl1^2+11520*m*xl2*xl1)*... x(i)^7+(17280*m*xl1^3-17280*m*xl2*xl1^2)*x(i)^6+(-11520*m*xl1^4+11520*m*xl2*... xl1^3+6*hbar^2)*x(i)^5+(2880*m*xl1^5-2880*m*xl2*xl1^4-60*hbar^2*xl1+54*hbar^2*... xl2)*x(i)^4+(108*hbar^2*xl1^2-108*hbar^2*xl2*xl1)*x(i)^3+(-132*hbar^2*xl1^3+132*... hbar^2*xl2*xl1^2)*x(i)^2+(78*hbar^2*xl1^4-78*hbar^2*xl2*xl1^3)*x(i)+(-18*hbar^2*... xl1^5+18*hbar^2*xl2*xl1^4)))*xr2+(((-2880*m*xl1+2880*m*xl2)*x(i)^4+(11520*m*... xl1^2-11520*m*xl2*xl1)*x(i)^3+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^2+(11520*... m*xl1^4-11520*m*xl2*xl1^3-6*hbar^2)*x(i)+(-2880*m*xl1^5+2880*m*xl2*xl1^4+24*... hbar^2*xl1-18*hbar^2*xl2))*xr1^5+((11520*m*xl1-11520*m*xl2)*x(i)^5+(-46080*m*... xl1^2+46080*m*xl2*xl1)*x(i)^4+(69120*m*xl1^3-69120*m*xl2*xl1^2)*x(i)^3+(-46080*... m*xl1^4+46080*m*xl2*xl1^3+24*hbar^2)*x(i)^2+(11520*m*xl1^5-11520*m*xl2*... xl1^4-102*hbar^2*xl1+78*hbar^2*xl2)*x(i)+(6*hbar^2*xl1^2-6*hbar^2*xl2*xl1))*... xr1^4+((-17280*m*xl1+17280*m*xl2)*x(i)^6+(69120*m*xl1^2-69120*m*xl2*xl1)*... x(i)^5+(-103680*m*xl1^3+103680*m*xl2*xl1^2)*x(i)^4+(69120*m*xl1^4-69120*m*xl2*... xl1^3-36*hbar^2)*x(i)^3+(-17280*m*xl1^5+17280*m*xl2*xl1^4+168*hbar^2*xl1-132*... hbar^2*xl2)*x(i)^2+(-30*hbar^2*xl1^2+30*hbar^2*xl2*xl1)*x(i)+(6*hbar^2*xl1^3-6*... hbar^2*xl2*xl1^2))*xr1^3+((11520*m*xl1-11520*m*xl2)*x(i)^7+(-46080*m*... xl1^2+46080*m*xl2*xl1)*x(i)^6+(69120*m*xl1^3-69120*m*xl2*xl1^2)*x(i)^5+(-46080*... m*xl1^4+46080*m*xl2*xl1^3+24*hbar^2)*x(i)^4+(11520*m*xl1^5-11520*m*xl2*... xl1^4-132*hbar^2*xl1+108*hbar^2*xl2)*x(i)^3+(60*hbar^2*xl1^2-60*hbar^2*xl2*xl1)*... x(i)^2+(-30*hbar^2*xl1^3+30*hbar^2*xl2*xl1^2)*x(i)+(6*hbar^2*xl1^4-6*hbar^2*xl2*... xl1^3))*xr1^2+((-2880*m*xl1+2880*m*xl2)*x(i)^8+(11520*m*xl1^2-11520*m*xl2*xl1)*... x(i)^7+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^6+(11520*m*xl1^4-11520*m*xl2*... xl1^3-6*hbar^2)*x(i)^5+(-2880*m*xl1^5+2880*m*xl2*xl1^4+66*hbar^2*xl1-60*hbar^2*... xl2)*x(i)^4+(-132*hbar^2*xl1^2+132*hbar^2*xl2*xl1)*x(i)^3+(168*hbar^2*xl1^3-168*... hbar^2*xl2*xl1^2)*x(i)^2+(-102*hbar^2*xl1^4+102*hbar^2*xl2*xl1^3)*x(i)+(24*... hbar^2*xl1^5-24*hbar^2*xl2*xl1^4))*xr1+((-6*hbar^2*xl1+6*hbar^2*xl2)*x(i)^5+(24*... hbar^2*xl1^2-24*hbar^2*xl2*xl1)*x(i)^4+(-36*hbar^2*xl1^3+36*hbar^2*xl2*xl1^2)*... x(i)^3+(24*hbar^2*xl1^4-24*hbar^2*xl2*xl1^3)*x(i)^2+(-6*hbar^2*xl1^5+6*hbar^2*... xl2*xl1^4)*x(i))))*exp(5*x(i)-5*c)^4+((((-720*m*xl1+720*m*xl2)*x(i)^4+(2880*m*... xl1^2-2880*m*xl2*xl1)*x(i)^3+(-4320*m*xl1^3+4320*m*xl2*xl1^2)*x(i)^2+(2880*m*... xl1^4-2880*m*xl2*xl1^3+4*hbar^2)*x(i)+(-720*m*xl1^5+720*m*xl2*xl1^4-16*hbar^2*... xl1+12*hbar^2*xl2))*xr1^4+((2880*m*xl1-2880*m*xl2)*x(i)^5+(-11520*m*xl1^2+11520*... m*xl2*xl1)*x(i)^4+(17280*m*xl1^3-17280*m*xl2*xl1^2)*x(i)^3+(-11520*m*... xl1^4+11520*m*xl2*xl1^3-16*hbar^2)*x(i)^2+(2880*m*xl1^5-2880*m*xl2*xl1^4+68*... hbar^2*xl1-52*hbar^2*xl2)*x(i)+(-4*hbar^2*xl1^2+4*hbar^2*xl2*xl1))*... xr1^3+((-4320*m*xl1+4320*m*xl2)*x(i)^6+(17280*m*xl1^2-17280*m*xl2*xl1)*... x(i)^5+(-25920*m*xl1^3+25920*m*xl2*xl1^2)*x(i)^4+(17280*m*xl1^4-17280*m*xl2*... xl1^3+24*hbar^2)*x(i)^3+(-4320*m*xl1^5+4320*m*xl2*xl1^4-112*hbar^2*xl1+88*... hbar^2*xl2)*x(i)^2+(20*hbar^2*xl1^2-20*hbar^2*xl2*xl1)*x(i)+(-4*hbar^2*xl1^3+4*... hbar^2*xl2*xl1^2))*xr1^2+((2880*m*xl1-2880*m*xl2)*x(i)^7+(-11520*m*xl1^2+11520*... m*xl2*xl1)*x(i)^6+(17280*m*xl1^3-17280*m*xl2*xl1^2)*x(i)^5+(-11520*m*... xl1^4+11520*m*xl2*xl1^3-16*hbar^2)*x(i)^4+(2880*m*xl1^5-2880*m*xl2*xl1^4+88*... hbar^2*xl1-72*hbar^2*xl2)*x(i)^3+(-40*hbar^2*xl1^2+40*hbar^2*xl2*xl1)*... x(i)^2+(20*hbar^2*xl1^3-20*hbar^2*xl2*xl1^2)*x(i)+(-4*hbar^2*xl1^4+4*hbar^2*xl2*... xl1^3))*xr1+((-720*m*xl1+720*m*xl2)*x(i)^8+(2880*m*xl1^2-2880*m*xl2*xl1)*... x(i)^7+(-4320*m*xl1^3+4320*m*xl2*xl1^2)*x(i)^6+(2880*m*xl1^4-2880*m*xl2*xl1^3+4*... hbar^2)*x(i)^5+(-720*m*xl1^5+720*m*xl2*xl1^4-40*hbar^2*xl1+36*hbar^2*xl2)*... x(i)^4+(72*hbar^2*xl1^2-72*hbar^2*xl2*xl1)*x(i)^3+(-88*hbar^2*xl1^3+88*hbar^2*... xl2*xl1^2)*x(i)^2+(52*hbar^2*xl1^4-52*hbar^2*xl2*xl1^3)*x(i)+(-12*hbar^2*... xl1^5+12*hbar^2*xl2*xl1^4)))*xr2+(((720*m*xl1-720*m*xl2)*x(i)^4+(-2880*m*... xl1^2+2880*m*xl2*xl1)*x(i)^3+(4320*m*xl1^3-4320*m*xl2*xl1^2)*x(i)^2+(-2880*m*... xl1^4+2880*m*xl2*xl1^3-4*hbar^2)*x(i)+(720*m*xl1^5-720*m*xl2*xl1^4+16*hbar^2*... xl1-12*hbar^2*xl2))*xr1^5+((-2880*m*xl1+2880*m*xl2)*x(i)^5+(11520*m*xl1^2-11520*... m*xl2*xl1)*x(i)^4+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^3+(11520*m*... xl1^4-11520*m*xl2*xl1^3+16*hbar^2)*x(i)^2+(-2880*m*xl1^5+2880*m*xl2*xl1^4-68*... hbar^2*xl1+52*hbar^2*xl2)*x(i)+(4*hbar^2*xl1^2-4*hbar^2*xl2*xl1))*xr1^4+((4320*... m*xl1-4320*m*xl2)*x(i)^6+(-17280*m*xl1^2+17280*m*xl2*xl1)*x(i)^5+(25920*m*... xl1^3-25920*m*xl2*xl1^2)*x(i)^4+(-17280*m*xl1^4+17280*m*xl2*xl1^3-24*hbar^2)*... x(i)^3+(4320*m*xl1^5-4320*m*xl2*xl1^4+112*hbar^2*xl1-88*hbar^2*xl2)*x(i)^2+(-20*... hbar^2*xl1^2+20*hbar^2*xl2*xl1)*x(i)+(4*hbar^2*xl1^3-4*hbar^2*xl2*xl1^2))*... xr1^3+((-2880*m*xl1+2880*m*xl2)*x(i)^7+(11520*m*xl1^2-11520*m*xl2*xl1)*... x(i)^6+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^5+(11520*m*xl1^4-11520*m*xl2*... xl1^3+16*hbar^2)*x(i)^4+(-2880*m*xl1^5+2880*m*xl2*xl1^4-88*hbar^2*xl1+72*hbar^2*... xl2)*x(i)^3+(40*hbar^2*xl1^2-40*hbar^2*xl2*xl1)*x(i)^2+(-20*hbar^2*xl1^3+20*... hbar^2*xl2*xl1^2)*x(i)+(4*hbar^2*xl1^4-4*hbar^2*xl2*xl1^3))*xr1^2+((720*m*... xl1-720*m*xl2)*x(i)^8+(-2880*m*xl1^2+2880*m*xl2*xl1)*x(i)^7+(4320*m*xl1^3-4320*... m*xl2*xl1^2)*x(i)^6+(-2880*m*xl1^4+2880*m*xl2*xl1^3-4*hbar^2)*x(i)^5+(720*m*... xl1^5-720*m*xl2*xl1^4+44*hbar^2*xl1-40*hbar^2*xl2)*x(i)^4+(-88*hbar^2*xl1^2+88*... hbar^2*xl2*xl1)*x(i)^3+(112*hbar^2*xl1^3-112*hbar^2*xl2*xl1^2)*x(i)^2+(-68*... hbar^2*xl1^4+68*hbar^2*xl2*xl1^3)*x(i)+(16*hbar^2*xl1^5-16*hbar^2*xl2*xl1^4))*... xr1+((-4*hbar^2*xl1+4*hbar^2*xl2)*x(i)^5+(16*hbar^2*xl1^2-16*hbar^2*xl2*xl1)*... x(i)^4+(-24*hbar^2*xl1^3+24*hbar^2*xl2*xl1^2)*x(i)^3+(16*hbar^2*xl1^4-16*hbar^2*... xl2*xl1^3)*x(i)^2+(-4*hbar^2*xl1^5+4*hbar^2*xl2*xl1^4)*x(i))))*exp(5*x(i)-5*... c)^2+(((hbar^2*x(i)+(-4*hbar^2*xl1+3*hbar^2*xl2))*xr1^4+(-4*hbar^2*x(i)^2+(17*... hbar^2*xl1-13*hbar^2*xl2)*x(i)+(-1*hbar^2*xl1^2+hbar^2*xl2*xl1))*xr1^3+(6*... hbar^2*x(i)^3+(-28*hbar^2*xl1+22*hbar^2*xl2)*x(i)^2+(5*hbar^2*xl1^2-5*hbar^2*... xl2*xl1)*x(i)+(-1*hbar^2*xl1^3+hbar^2*xl2*xl1^2))*xr1^2+(-4*hbar^2*x(i)^4+(22*... hbar^2*xl1-18*hbar^2*xl2)*x(i)^3+(-10*hbar^2*xl1^2+10*hbar^2*xl2*xl1)*x(i)^2+(5*... hbar^2*xl1^3-5*hbar^2*xl2*xl1^2)*x(i)+(-1*hbar^2*xl1^4+hbar^2*xl2*xl1^3))*... xr1+(hbar^2*x(i)^5+(-10*hbar^2*xl1+9*hbar^2*xl2)*x(i)^4+(18*hbar^2*xl1^2-18*... hbar^2*xl2*xl1)*x(i)^3+(-22*hbar^2*xl1^3+22*hbar^2*xl2*xl1^2)*x(i)^2+(13*hbar^2*... xl1^4-13*hbar^2*xl2*xl1^3)*x(i)+(-3*hbar^2*xl1^5+3*hbar^2*xl2*xl1^4)))*xr2+((-1*... hbar^2*x(i)+(4*hbar^2*xl1-3*hbar^2*xl2))*xr1^5+(4*hbar^2*x(i)^2+(-17*hbar^2*... xl1+13*hbar^2*xl2)*x(i)+(hbar^2*xl1^2-1*hbar^2*xl2*xl1))*xr1^4+(-6*hbar^2*... x(i)^3+(28*hbar^2*xl1-22*hbar^2*xl2)*x(i)^2+(-5*hbar^2*xl1^2+5*hbar^2*xl2*xl1)*... x(i)+(hbar^2*xl1^3-1*hbar^2*xl2*xl1^2))*xr1^3+(4*hbar^2*x(i)^4+(-22*hbar^2*... xl1+18*hbar^2*xl2)*x(i)^3+(10*hbar^2*xl1^2-10*hbar^2*xl2*xl1)*x(i)^2+(-5*hbar^2*... xl1^3+5*hbar^2*xl2*xl1^2)*x(i)+(hbar^2*xl1^4-1*hbar^2*xl2*xl1^3))*xr1^2+(-1*... hbar^2*x(i)^5+(11*hbar^2*xl1-10*hbar^2*xl2)*x(i)^4+(-22*hbar^2*xl1^2+22*hbar^2*... xl2*xl1)*x(i)^3+(28*hbar^2*xl1^3-28*hbar^2*xl2*xl1^2)*x(i)^2+(-17*hbar^2*... xl1^4+17*hbar^2*xl2*xl1^3)*x(i)+(4*hbar^2*xl1^5-4*hbar^2*xl2*xl1^4))*xr1+((-1*... hbar^2*xl1+hbar^2*xl2)*x(i)^5+(4*hbar^2*xl1^2-4*hbar^2*xl2*xl1)*x(i)^4+(-6*... hbar^2*xl1^3+6*hbar^2*xl2*xl1^2)*x(i)^3+(4*hbar^2*xl1^4-4*hbar^2*xl2*xl1^3)*... x(i)^2+(-1*hbar^2*xl1^5+hbar^2*xl2*xl1^4)*x(i)))))/(((((2*m^2*xl1-2*m^2*xl2)*... x(i)^4+(-8*m^2*xl1^2+8*m^2*xl2*xl1)*x(i)^3+(12*m^2*xl1^3-12*m^2*xl2*xl1^2)*... x(i)^2+(-8*m^2*xl1^4+8*m^2*xl2*xl1^3)*x(i)+(2*m^2*xl1^5-2*m^2*xl2*xl1^4))*... xr1^4+((-8*m^2*xl1+8*m^2*xl2)*x(i)^5+(32*m^2*xl1^2-32*m^2*xl2*xl1)*x(i)^4+(-48*... m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^3+(32*m^2*xl1^4-32*m^2*xl2*xl1^3)*x(i)^2+(-8*... m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i))*xr1^3+((12*m^2*xl1-12*m^2*xl2)*x(i)^6+(-48*m^2*... xl1^2+48*m^2*xl2*xl1)*x(i)^5+(72*m^2*xl1^3-72*m^2*xl2*xl1^2)*x(i)^4+(-48*m^2*... xl1^4+48*m^2*xl2*xl1^3)*x(i)^3+(12*m^2*xl1^5-12*m^2*xl2*xl1^4)*x(i)^2)*... xr1^2+((-8*m^2*xl1+8*m^2*xl2)*x(i)^7+(32*m^2*xl1^2-32*m^2*xl2*xl1)*x(i)^6+(-48*... m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^5+(32*m^2*xl1^4-32*m^2*xl2*xl1^3)*x(i)^4+(-8*... m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i)^3)*xr1+((2*m^2*xl1-2*m^2*xl2)*x(i)^8+(-8*m^2*... xl1^2+8*m^2*xl2*xl1)*x(i)^7+(12*m^2*xl1^3-12*m^2*xl2*xl1^2)*x(i)^6+(-8*m^2*... xl1^4+8*m^2*xl2*xl1^3)*x(i)^5+(2*m^2*xl1^5-2*m^2*xl2*xl1^4)*x(i)^4))*xr2+(((-2*... m^2*xl1+2*m^2*xl2)*x(i)^4+(8*m^2*xl1^2-8*m^2*xl2*xl1)*x(i)^3+(-12*m^2*xl1^3+12*... m^2*xl2*xl1^2)*x(i)^2+(8*m^2*xl1^4-8*m^2*xl2*xl1^3)*x(i)+(-2*m^2*xl1^5+2*m^2*... xl2*xl1^4))*xr1^5+((8*m^2*xl1-8*m^2*xl2)*x(i)^5+(-32*m^2*xl1^2+32*m^2*xl2*xl1)*... x(i)^4+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^3+(-32*m^2*xl1^4+32*m^2*xl2*xl1^3)*... x(i)^2+(8*m^2*xl1^5-8*m^2*xl2*xl1^4)*x(i))*xr1^4+((-12*m^2*xl1+12*m^2*xl2)*... x(i)^6+(48*m^2*xl1^2-48*m^2*xl2*xl1)*x(i)^5+(-72*m^2*xl1^3+72*m^2*xl2*xl1^2)*... x(i)^4+(48*m^2*xl1^4-48*m^2*xl2*xl1^3)*x(i)^3+(-12*m^2*xl1^5+12*m^2*xl2*xl1^4)*... x(i)^2)*xr1^3+((8*m^2*xl1-8*m^2*xl2)*x(i)^7+(-32*m^2*xl1^2+32*m^2*xl2*xl1)*... x(i)^6+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^5+(-32*m^2*xl1^4+32*m^2*xl2*xl1^3)*... x(i)^4+(8*m^2*xl1^5-8*m^2*xl2*xl1^4)*x(i)^3)*xr1^2+((-2*m^2*xl1+2*m^2*xl2)*... x(i)^8+(8*m^2*xl1^2-8*m^2*xl2*xl1)*x(i)^7+(-12*m^2*xl1^3+12*m^2*xl2*xl1^2)*... x(i)^6+(8*m^2*xl1^4-8*m^2*xl2*xl1^3)*x(i)^5+(-2*m^2*xl1^5+2*m^2*xl2*xl1^4)*... x(i)^4)*xr1))*exp(5*x(i)-5*c)^8+((((8*m^2*xl1-8*m^2*xl2)*x(i)^4+(-32*m^2*... xl1^2+32*m^2*xl2*xl1)*x(i)^3+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^2+(-32*m^2*... xl1^4+32*m^2*xl2*xl1^3)*x(i)+(8*m^2*xl1^5-8*m^2*xl2*xl1^4))*xr1^4+((-32*m^2*... xl1+32*m^2*xl2)*x(i)^5+(128*m^2*xl1^2-128*m^2*xl2*xl1)*x(i)^4+(-192*m^2*... xl1^3+192*m^2*xl2*xl1^2)*x(i)^3+(128*m^2*xl1^4-128*m^2*xl2*xl1^3)*x(i)^2+(-32*... m^2*xl1^5+32*m^2*xl2*xl1^4)*x(i))*xr1^3+((48*m^2*xl1-48*m^2*xl2)*x(i)^6+(-192*... m^2*xl1^2+192*m^2*xl2*xl1)*x(i)^5+(288*m^2*xl1^3-288*m^2*xl2*xl1^2)*... x(i)^4+(-192*m^2*xl1^4+192*m^2*xl2*xl1^3)*x(i)^3+(48*m^2*xl1^5-48*m^2*xl2*... xl1^4)*x(i)^2)*xr1^2+((-32*m^2*xl1+32*m^2*xl2)*x(i)^7+(128*m^2*xl1^2-128*m^2*... xl2*xl1)*x(i)^6+(-192*m^2*xl1^3+192*m^2*xl2*xl1^2)*x(i)^5+(128*m^2*xl1^4-128*... m^2*xl2*xl1^3)*x(i)^4+(-32*m^2*xl1^5+32*m^2*xl2*xl1^4)*x(i)^3)*xr1+((8*m^2*... xl1-8*m^2*xl2)*x(i)^8+(-32*m^2*xl1^2+32*m^2*xl2*xl1)*x(i)^7+(48*m^2*xl1^3-48*... m^2*xl2*xl1^2)*x(i)^6+(-32*m^2*xl1^4+32*m^2*xl2*xl1^3)*x(i)^5+(8*m^2*xl1^5-8*... m^2*xl2*xl1^4)*x(i)^4))*xr2+(((-8*m^2*xl1+8*m^2*xl2)*x(i)^4+(32*m^2*xl1^2-32*... m^2*xl2*xl1)*x(i)^3+(-48*m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^2+(32*m^2*xl1^4-32*... m^2*xl2*xl1^3)*x(i)+(-8*m^2*xl1^5+8*m^2*xl2*xl1^4))*xr1^5+((32*m^2*xl1-32*m^2*... xl2)*x(i)^5+(-128*m^2*xl1^2+128*m^2*xl2*xl1)*x(i)^4+(192*m^2*xl1^3-192*m^2*xl2*... xl1^2)*x(i)^3+(-128*m^2*xl1^4+128*m^2*xl2*xl1^3)*x(i)^2+(32*m^2*xl1^5-32*m^2*... xl2*xl1^4)*x(i))*xr1^4+((-48*m^2*xl1+48*m^2*xl2)*x(i)^6+(192*m^2*xl1^2-192*m^2*... xl2*xl1)*x(i)^5+(-288*m^2*xl1^3+288*m^2*xl2*xl1^2)*x(i)^4+(192*m^2*xl1^4-192*... m^2*xl2*xl1^3)*x(i)^3+(-48*m^2*xl1^5+48*m^2*xl2*xl1^4)*x(i)^2)*xr1^3+((32*m^2*... xl1-32*m^2*xl2)*x(i)^7+(-128*m^2*xl1^2+128*m^2*xl2*xl1)*x(i)^6+(192*m^2*... xl1^3-192*m^2*xl2*xl1^2)*x(i)^5+(-128*m^2*xl1^4+128*m^2*xl2*xl1^3)*x(i)^4+(32*... m^2*xl1^5-32*m^2*xl2*xl1^4)*x(i)^3)*xr1^2+((-8*m^2*xl1+8*m^2*xl2)*x(i)^8+(32*... m^2*xl1^2-32*m^2*xl2*xl1)*x(i)^7+(-48*m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^6+(32*... m^2*xl1^4-32*m^2*xl2*xl1^3)*x(i)^5+(-8*m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i)^4)*xr1))*... exp(5*x(i)-5*c)^6+((((12*m^2*xl1-12*m^2*xl2)*x(i)^4+(-48*m^2*xl1^2+48*m^2*xl2*... xl1)*x(i)^3+(72*m^2*xl1^3-72*m^2*xl2*xl1^2)*x(i)^2+(-48*m^2*xl1^4+48*m^2*xl2*... xl1^3)*x(i)+(12*m^2*xl1^5-12*m^2*xl2*xl1^4))*xr1^4+((-48*m^2*xl1+48*m^2*xl2)*... x(i)^5+(192*m^2*xl1^2-192*m^2*xl2*xl1)*x(i)^4+(-288*m^2*xl1^3+288*m^2*xl2*... xl1^2)*x(i)^3+(192*m^2*xl1^4-192*m^2*xl2*xl1^3)*x(i)^2+(-48*m^2*xl1^5+48*m^2*... xl2*xl1^4)*x(i))*xr1^3+((72*m^2*xl1-72*m^2*xl2)*x(i)^6+(-288*m^2*xl1^2+288*m^2*... xl2*xl1)*x(i)^5+(432*m^2*xl1^3-432*m^2*xl2*xl1^2)*x(i)^4+(-288*m^2*xl1^4+288*... m^2*xl2*xl1^3)*x(i)^3+(72*m^2*xl1^5-72*m^2*xl2*xl1^4)*x(i)^2)*xr1^2+((-48*m^2*... xl1+48*m^2*xl2)*x(i)^7+(192*m^2*xl1^2-192*m^2*xl2*xl1)*x(i)^6+(-288*m^2*... xl1^3+288*m^2*xl2*xl1^2)*x(i)^5+(192*m^2*xl1^4-192*m^2*xl2*xl1^3)*x(i)^4+(-48*... m^2*xl1^5+48*m^2*xl2*xl1^4)*x(i)^3)*xr1+((12*m^2*xl1-12*m^2*xl2)*x(i)^8+(-48*... m^2*xl1^2+48*m^2*xl2*xl1)*x(i)^7+(72*m^2*xl1^3-72*m^2*xl2*xl1^2)*x(i)^6+(-48*... m^2*xl1^4+48*m^2*xl2*xl1^3)*x(i)^5+(12*m^2*xl1^5-12*m^2*xl2*xl1^4)*x(i)^4))*... xr2+(((-12*m^2*xl1+12*m^2*xl2)*x(i)^4+(48*m^2*xl1^2-48*m^2*xl2*xl1)*x(i)^3+(-72*... m^2*xl1^3+72*m^2*xl2*xl1^2)*x(i)^2+(48*m^2*xl1^4-48*m^2*xl2*xl1^3)*x(i)+(-12*... m^2*xl1^5+12*m^2*xl2*xl1^4))*xr1^5+((48*m^2*xl1-48*m^2*xl2)*x(i)^5+(-192*m^2*... xl1^2+192*m^2*xl2*xl1)*x(i)^4+(288*m^2*xl1^3-288*m^2*xl2*xl1^2)*x(i)^3+(-192*... m^2*xl1^4+192*m^2*xl2*xl1^3)*x(i)^2+(48*m^2*xl1^5-48*m^2*xl2*xl1^4)*x(i))*... xr1^4+((-72*m^2*xl1+72*m^2*xl2)*x(i)^6+(288*m^2*xl1^2-288*m^2*xl2*xl1)*... x(i)^5+(-432*m^2*xl1^3+432*m^2*xl2*xl1^2)*x(i)^4+(288*m^2*xl1^4-288*m^2*xl2*... xl1^3)*x(i)^3+(-72*m^2*xl1^5+72*m^2*xl2*xl1^4)*x(i)^2)*xr1^3+((48*m^2*xl1-48*... m^2*xl2)*x(i)^7+(-192*m^2*xl1^2+192*m^2*xl2*xl1)*x(i)^6+(288*m^2*xl1^3-288*m^2*... xl2*xl1^2)*x(i)^5+(-192*m^2*xl1^4+192*m^2*xl2*xl1^3)*x(i)^4+(48*m^2*xl1^5-48*... m^2*xl2*xl1^4)*x(i)^3)*xr1^2+((-12*m^2*xl1+12*m^2*xl2)*x(i)^8+(48*m^2*xl1^2-48*... m^2*xl2*xl1)*x(i)^7+(-72*m^2*xl1^3+72*m^2*xl2*xl1^2)*x(i)^6+(48*m^2*xl1^4-48*... m^2*xl2*xl1^3)*x(i)^5+(-12*m^2*xl1^5+12*m^2*xl2*xl1^4)*x(i)^4)*xr1))*exp(5*... x(i)-5*c)^4+((((8*m^2*xl1-8*m^2*xl2)*x(i)^4+(-32*m^2*xl1^2+32*m^2*xl2*xl1)*... x(i)^3+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^2+(-32*m^2*xl1^4+32*m^2*xl2*xl1^3)*... x(i)+(8*m^2*xl1^5-8*m^2*xl2*xl1^4))*xr1^4+((-32*m^2*xl1+32*m^2*xl2)*x(i)^5+(128*... m^2*xl1^2-128*m^2*xl2*xl1)*x(i)^4+(-192*m^2*xl1^3+192*m^2*xl2*xl1^2)*... x(i)^3+(128*m^2*xl1^4-128*m^2*xl2*xl1^3)*x(i)^2+(-32*m^2*xl1^5+32*m^2*xl2*... xl1^4)*x(i))*xr1^3+((48*m^2*xl1-48*m^2*xl2)*x(i)^6+(-192*m^2*xl1^2+192*m^2*xl2*... xl1)*x(i)^5+(288*m^2*xl1^3-288*m^2*xl2*xl1^2)*x(i)^4+(-192*m^2*xl1^4+192*m^2*... xl2*xl1^3)*x(i)^3+(48*m^2*xl1^5-48*m^2*xl2*xl1^4)*x(i)^2)*xr1^2+((-32*m^2*... xl1+32*m^2*xl2)*x(i)^7+(128*m^2*xl1^2-128*m^2*xl2*xl1)*x(i)^6+(-192*m^2*... xl1^3+192*m^2*xl2*xl1^2)*x(i)^5+(128*m^2*xl1^4-128*m^2*xl2*xl1^3)*x(i)^4+(-32*... m^2*xl1^5+32*m^2*xl2*xl1^4)*x(i)^3)*xr1+((8*m^2*xl1-8*m^2*xl2)*x(i)^8+(-32*m^2*... xl1^2+32*m^2*xl2*xl1)*x(i)^7+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^6+(-32*m^2*... xl1^4+32*m^2*xl2*xl1^3)*x(i)^5+(8*m^2*xl1^5-8*m^2*xl2*xl1^4)*x(i)^4))*xr2+(((-8*... m^2*xl1+8*m^2*xl2)*x(i)^4+(32*m^2*xl1^2-32*m^2*xl2*xl1)*x(i)^3+(-48*m^2*... xl1^3+48*m^2*xl2*xl1^2)*x(i)^2+(32*m^2*xl1^4-32*m^2*xl2*xl1^3)*x(i)+(-8*m^2*... xl1^5+8*m^2*xl2*xl1^4))*xr1^5+((32*m^2*xl1-32*m^2*xl2)*x(i)^5+(-128*m^2*... xl1^2+128*m^2*xl2*xl1)*x(i)^4+(192*m^2*xl1^3-192*m^2*xl2*xl1^2)*x(i)^3+(-128*... m^2*xl1^4+128*m^2*xl2*xl1^3)*x(i)^2+(32*m^2*xl1^5-32*m^2*xl2*xl1^4)*x(i))*... xr1^4+((-48*m^2*xl1+48*m^2*xl2)*x(i)^6+(192*m^2*xl1^2-192*m^2*xl2*xl1)*... x(i)^5+(-288*m^2*xl1^3+288*m^2*xl2*xl1^2)*x(i)^4+(192*m^2*xl1^4-192*m^2*xl2*... xl1^3)*x(i)^3+(-48*m^2*xl1^5+48*m^2*xl2*xl1^4)*x(i)^2)*xr1^3+((32*m^2*xl1-32*... m^2*xl2)*x(i)^7+(-128*m^2*xl1^2+128*m^2*xl2*xl1)*x(i)^6+(192*m^2*xl1^3-192*m^2*... xl2*xl1^2)*x(i)^5+(-128*m^2*xl1^4+128*m^2*xl2*xl1^3)*x(i)^4+(32*m^2*xl1^5-32*... m^2*xl2*xl1^4)*x(i)^3)*xr1^2+((-8*m^2*xl1+8*m^2*xl2)*x(i)^8+(32*m^2*xl1^2-32*... m^2*xl2*xl1)*x(i)^7+(-48*m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^6+(32*m^2*xl1^4-32*... m^2*xl2*xl1^3)*x(i)^5+(-8*m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i)^4)*xr1))*exp(5*x(i)-5*... c)^2+((((2*m^2*xl1-2*m^2*xl2)*x(i)^4+(-8*m^2*xl1^2+8*m^2*xl2*xl1)*x(i)^3+(12*... m^2*xl1^3-12*m^2*xl2*xl1^2)*x(i)^2+(-8*m^2*xl1^4+8*m^2*xl2*xl1^3)*x(i)+(2*m^2*... xl1^5-2*m^2*xl2*xl1^4))*xr1^4+((-8*m^2*xl1+8*m^2*xl2)*x(i)^5+(32*m^2*xl1^2-32*... m^2*xl2*xl1)*x(i)^4+(-48*m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^3+(32*m^2*xl1^4-32*... m^2*xl2*xl1^3)*x(i)^2+(-8*m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i))*xr1^3+((12*m^2*... xl1-12*m^2*xl2)*x(i)^6+(-48*m^2*xl1^2+48*m^2*xl2*xl1)*x(i)^5+(72*m^2*xl1^3-72*... m^2*xl2*xl1^2)*x(i)^4+(-48*m^2*xl1^4+48*m^2*xl2*xl1^3)*x(i)^3+(12*m^2*xl1^5-12*... m^2*xl2*xl1^4)*x(i)^2)*xr1^2+((-8*m^2*xl1+8*m^2*xl2)*x(i)^7+(32*m^2*xl1^2-32*... m^2*xl2*xl1)*x(i)^6+(-48*m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^5+(32*m^2*xl1^4-32*... m^2*xl2*xl1^3)*x(i)^4+(-8*m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i)^3)*xr1+((2*m^2*xl1-2*... m^2*xl2)*x(i)^8+(-8*m^2*xl1^2+8*m^2*xl2*xl1)*x(i)^7+(12*m^2*xl1^3-12*m^2*xl2*... xl1^2)*x(i)^6+(-8*m^2*xl1^4+8*m^2*xl2*xl1^3)*x(i)^5+(2*m^2*xl1^5-2*m^2*xl2*... xl1^4)*x(i)^4))*xr2+(((-2*m^2*xl1+2*m^2*xl2)*x(i)^4+(8*m^2*xl1^2-8*m^2*xl2*xl1)*... x(i)^3+(-12*m^2*xl1^3+12*m^2*xl2*xl1^2)*x(i)^2+(8*m^2*xl1^4-8*m^2*xl2*xl1^3)*... x(i)+(-2*m^2*xl1^5+2*m^2*xl2*xl1^4))*xr1^5+((8*m^2*xl1-8*m^2*xl2)*x(i)^5+(-32*... m^2*xl1^2+32*m^2*xl2*xl1)*x(i)^4+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^3+(-32*... m^2*xl1^4+32*m^2*xl2*xl1^3)*x(i)^2+(8*m^2*xl1^5-8*m^2*xl2*xl1^4)*x(i))*... xr1^4+((-12*m^2*xl1+12*m^2*xl2)*x(i)^6+(48*m^2*xl1^2-48*m^2*xl2*xl1)*... x(i)^5+(-72*m^2*xl1^3+72*m^2*xl2*xl1^2)*x(i)^4+(48*m^2*xl1^4-48*m^2*xl2*xl1^3)*... x(i)^3+(-12*m^2*xl1^5+12*m^2*xl2*xl1^4)*x(i)^2)*xr1^3+((8*m^2*xl1-8*m^2*xl2)*... x(i)^7+(-32*m^2*xl1^2+32*m^2*xl2*xl1)*x(i)^6+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*... x(i)^5+(-32*m^2*xl1^4+32*m^2*xl2*xl1^3)*x(i)^4+(8*m^2*xl1^5-8*m^2*xl2*xl1^4)*... x(i)^3)*xr1^2+((-2*m^2*xl1+2*m^2*xl2)*x(i)^8+(8*m^2*xl1^2-8*m^2*xl2*xl1)*... x(i)^7+(-12*m^2*xl1^3+12*m^2*xl2*xl1^2)*x(i)^6+(8*m^2*xl1^4-8*m^2*xl2*xl1^3)*... x(i)^5+(-2*m^2*xl1^5+2*m^2*xl2*xl1^4)*x(i)^4)*xr1))); endfunction
%sage print fricas2octave('Jp1','jac1(i,1)')

function result=Jp1(x,i) global l2 l1 p1 p2; global N m hbar a b c; if (i-2<1) xl2=l2; else xl2=x(i-2); endif; if (i-1<1) xl1=l1; else xl1=x(i-1); endif; if (i+2>N) xr2=p2; else xr2=x(i+2); endif; if (i+1>N) xr1=p1; else xr1=x(i+1); endif; result=((hbar^2*xr1^2+(-4*hbar^2*x(i)+2*hbar^2*xl1)*xr1+(9*hbar^2*x(i)^2-14*hbar^2*xl1*... x(i)+6*hbar^2*xl1^2))*xr2^2+(-2*hbar^2*xr1^3+(8*hbar^2*x(i)-4*hbar^2*xl1)*... xr1^2+(-20*hbar^2*x(i)^2+32*hbar^2*xl1*x(i)-14*hbar^2*xl1^2)*xr1+(2*hbar^2*... x(i)^3-4*hbar^2*xl1*x(i)^2+2*hbar^2*xl1^2*x(i)))*xr2+(hbar^2*xr1^4+(-4*hbar^2*... x(i)+2*hbar^2*xl1)*xr1^3+(12*hbar^2*x(i)^2-20*hbar^2*xl1*x(i)+9*hbar^2*xl1^2)*... xr1^2+(-4*hbar^2*x(i)^3+8*hbar^2*xl1*x(i)^2-4*hbar^2*xl1^2*x(i))*xr1+(hbar^2*... x(i)^4-2*hbar^2*xl1*x(i)^3+hbar^2*xl1^2*x(i)^2)))/(((4*m^2*x(i)^2-8*m^2*xl1*... x(i)+4*m^2*xl1^2)*xr1^4+(-16*m^2*x(i)^3+32*m^2*xl1*x(i)^2-16*m^2*xl1^2*x(i))*... xr1^3+(24*m^2*x(i)^4-48*m^2*xl1*x(i)^3+24*m^2*xl1^2*x(i)^2)*xr1^2+(-16*m^2*... x(i)^5+32*m^2*xl1*x(i)^4-16*m^2*xl1^2*x(i)^3)*xr1+(4*m^2*x(i)^6-8*m^2*xl1*... x(i)^5+4*m^2*xl1^2*x(i)^4))*xr2^2+((-8*m^2*x(i)^2+16*m^2*xl1*x(i)-8*m^2*xl1^2)*... xr1^5+(32*m^2*x(i)^3-64*m^2*xl1*x(i)^2+32*m^2*xl1^2*x(i))*xr1^4+(-48*m^2*... x(i)^4+96*m^2*xl1*x(i)^3-48*m^2*xl1^2*x(i)^2)*xr1^3+(32*m^2*x(i)^5-64*m^2*xl1*... x(i)^4+32*m^2*xl1^2*x(i)^3)*xr1^2+(-8*m^2*x(i)^6+16*m^2*xl1*x(i)^5-8*m^2*xl1^2*... x(i)^4)*xr1)*xr2+((4*m^2*x(i)^2-8*m^2*xl1*x(i)+4*m^2*xl1^2)*xr1^6+(-16*m^2*... x(i)^3+32*m^2*xl1*x(i)^2-16*m^2*xl1^2*x(i))*xr1^5+(24*m^2*x(i)^4-48*m^2*xl1*... x(i)^3+24*m^2*xl1^2*x(i)^2)*xr1^4+(-16*m^2*x(i)^5+32*m^2*xl1*x(i)^4-16*m^2*... xl1^2*x(i)^3)*xr1^3+(4*m^2*x(i)^6-8*m^2*xl1*x(i)^5+4*m^2*xl1^2*x(i)^4)*xr1^2)); endfunction
%sage print fricas2octave('Jp2','jac1(i,2)')

function result=Jp2(x,i) global l2 l1 p1 p2; global N m hbar a b c; if (i-2<1) xl2=l2; else xl2=x(i-2); endif; if (i-1<1) xl1=l1; else xl1=x(i-1); endif; if (i+2>N) xr2=p2; else xr2=x(i+2); endif; if (i+1>N) xr1=p1; else xr1=x(i+1); endif; result=(-1*hbar^2)/((4*m^2*xr1^2-8*m^2*x(i)*xr1+4*m^2*x(i)^2)*xr2^2+(-8*m^2*xr1^3+16*... m^2*x(i)*xr1^2-8*m^2*x(i)^2*xr1)*xr2+(4*m^2*xr1^4-8*m^2*x(i)*xr1^3+4*m^2*x(i)^2*... xr1^2)); endfunction
function result=jac1(x,i,j)
global m hbar a;
if (j==-2)
result=Jm2(x,i);
elseif (j==-1)
result=Jm1(x,i);
elseif (j==0)
result=J0(x,i);
elseif (j==1)
result=Jp1(x,i);
elseif (j==2)
result=Jp2(x,i);
else
result=0;
endif
endfunction

jac1([X0;Y0],3,2)

ans = -2.50144
function result = jac0(x,i,j)
global N m hbar a;
if (i==j-N)
result = 1;
elseif (i<=N)
result = 0;
elseif (i<N+j-2)
result = 0;
else
result = jac1(x,i-N,i-N-j);
endif
endfunction

jac0([X0;Y0],15,5)

ans = 0
function jac=jacm(x,t)
global N
jac=zeros(2*N,2*N);
for i=1:2*N
for j=1:2*N
jac(i,j)=jac0(x,i,j);
endfor
endfor
endfunction


t0 = clock ();
trec = linspace (0, 10, 240)';

#Xrec = lsode (@deq, [X0;Y0], trec)';
Xrec = lsode ({@deq,@jacm}, [X0;Y0], trec)';
elapsed_time = etime (clock (), t0)
T=size(trec,1);
lsode_options()


elapsed_time = 108.378 Options for LSODE include: keyword value ------- ----- absolute tolerance 1.49012e-08 relative tolerance 1.49012e-08 integration method stiff initial step size -1 maximum order -1 maximum step size -1 minimum step size 0 step limit 100000
h = figure('visible', 'off');
plot(trec,Xrec(1:N,1:T));
saveas(h,"figure2.png");
%ylim([-2,2]);

%sage salvus.file("figure2.png") h = figure('visible', 'off');
%Plot the density at the end of the simulation
Dist=((Xrec(2:N,T)-Xrec(1:(N-1),T)));
Prob=(1/N)*Dist.^(-1);
plot(Xrec(2:N,T),Prob);
saveas(h,"figure3.png");
%xlim([-1.5,1.5]); ylim([0,1]);

%sage salvus.file("figure3.png")

h = figure('visible', 'off');
%Plot the density at 3/4 through the simulation
Tsamp=3*T/4;
Dist=((Xrec(2:N,Tsamp)-Xrec(1:(N-1),Tsamp)));
Prob=(1/N)*Dist.^(-1);

plot(Xrec(2:N,Tsamp),Prob);
saveas(h,"figure4.png");
%xlim([-1.5,1.5]); ylim([0,1]);

%sage salvus.file("figure4.png") h = figure('visible', 'off');
Dist=((Xrec(2:N,T)-Xrec(1:(N-1),T)));
Prob=(1/N)*Dist.^(-1);
plot(xvalues(2:N,1),Pin);
hold on
Dist=((Xrec(2:N,T)-Xrec(1:(N-1),T)));
Prob=(1/N)*Dist.^(-1);
plot(Xrec(2:N,T),Prob,"r");
saveas(h,"figure5.png");
hold off

%sage salvus.file("figure5.png") figure("visible","off");
moviename="MIW1Dv4.mp4";
movie("init",moviename)
for Tsamp=1:10:T
Dist=((Xrec(2:N,Tsamp)-Xrec(1:(N-1),Tsamp)));
Prob=Dist.^(-1);
plot(Xrec(2:N,Tsamp),Prob);
xlim([-20,20]); ylim([0,20]);
endfor
movie("close",moviename,4); %last arg is frame rate
T

T = 240
%sage salvus.file("MIW1Dv4.mp4")