CoCalc Shared FilesMIWsim1DjacobianV1b.sagewsOpen in CoCalc with one click!
Authors: David Cyganski, Bill Page
Views : 9
Description: Implementation of Wiseman Many-Interacting-Worlds simulation for 1D, 1 Particle systems

###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)=[x1(t)...xNJD(t)]Xp(t)=[x_1(t)...x_{N*J*D}(t)] is determined in general in MIW by system of NJDN*J*D coupled non-linear second order partial differential equations, where NN=number of worlds, JJ=number of particles, DD=number of dimensions. mkx¨kn(t)=fk(xn(t))+rkN(xn(t);Xp(t)) m_k*\ddot{x}_{kn}(t) = f_k(x_n(t)) + r_{kN}(x_n(t);Xp(t)) where kk indexes the particle and nn the world.

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

x˙kn(t)=ykn(t)y˙kn(t)=(1/mk)(fk(xn(t))+rkN(xn(t)) \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 2N2N 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=cx=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 ithi^{th} world: x(i)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

Perform the time updates

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]); movie("add",moviename); endfor movie("close",moviename,4); %last arg is frame rate T
T = 240
%sage salvus.file("MIW1Dv4.mp4")