## 5038 Reputation

9 years, 219 days

## this is the answer I meant to give...

but my original comment about not really understanding teh problem still stands.

However at least this answer is vaguely(?) appropriate for this question

 > restart: Severity:=[1,2,3]: Likelihood:=[1,2,3]: Impact:=LinearAlgebra:-OuterProductMatrix(Severity, Likelihood); Outcome:=Matrix( op(1, Impact),                  (i,j)-> if Impact[i,j] <= 3                          then "minor"                          elif Impact[i,j] <= 6                          then "moderate"                          else "serious"                          fi                );
 (1)
 >

## Maybe this will help...

I'm not really sure what you are trying to achieve.

The attached will

1. Plot five more or less random graphs
2. Extract the data from these plots and build them into a matrix where the first column contains x-values: subsequent five columns contain y-values for each of the five graphs in (1) above
3. Export this matrix (in csv format - many other formats are possible)
 > restart;   with(plots):   with(plottools): # # Define five more or less random functions to plot #   v:= Vector       ( 5,         [1.1*t, 1.2*t, 1.3*t, 1.4*t, 1.5*t]       ): # # Plot the five functions above #   g11:= plot         ( [ seq             ( v[j],               j=1..5             )           ],           t = 1 .. 2.9,           color = [red, green, blue, black, purple],           numpoints = 200         ); # # Extract the plot data from the above in the form of # a matrix. # # First column is x-values: subsequent columns are the # y-values associated with each of the individual plots #   M1:= Matrix        ( [ getdata(g11)[1][3][..,1],            seq            ( getdata(g11)[j][3][..,2],              j=1..5            )          ],          scan=rows        ); # # Export this matrix to a csv file on my desktop. # # 1. Obviously the OP will have to change the #    location of the output file to something #    appropriate for his/her machine/OS # 2. Output will be in plain vanilla 'csv' format. #    Many other output formats are possible. #   ExportMatrix   ( "C:/Users/TomLeslie/Desktop/plotData.csv",     M1,     target=csv   );
 (1)
 >
 >

## Undefined versus discontinuous...

You have to dra a distinction between whether a function is "undefined" or "discontinuous"

The simple Heaviside function, ie Heaviside(x) has the value 0, for x<0, and the value 1, for x>0. Note that

1. It is disconcinuous at x=0
2. It is undefined at x=0

However consider the "simple" piecewise function

f:=piecewise(x<0,0,
x=0, 1/2,
x>0, 1
)

1. It is discontinuous at x=0
2. It is defined for all x - including x=0. Note that the value at x=0 given in the above is "arbitrary". Other "quite sensible" values would be 0 or 1. In fact, even 100: this would still mean that the function is "defined" for all x - just "discontinuous" at x=0

In the attached I have shown what happens when you replace your "Heaviside" definitions with "piecewise" definitions.

Basically fsolve() never gives an answer, even although the piecewise function is everywhere defined.

On the other hand DirectSearch:SolveEquations() always gives an answer - althugh you may want to think about these carefully, given the reported value of the "residual"

 >
 > restart:
 > with(DirectSearch):
 > F:=x*(Heaviside(x-(S0))/(1-DD)+Heaviside(-x+S0)/(1-HH*DD));
 (1)
 > EQ:=F-E*Y;
 (2)
 > HH:=0.2; DD:=0.579350262599427; S0:=1e8; E:=70e9; Y:=0.001668993263529;
 (3)
 > EQ;
 (4)
 > plot(EQ, x=0..10^9, discont=true);
 > fsolve(EQ)
 (5)
 > SolveEquations(EQ)
 (6)
 > eval(EQ,x = 1.00000000000000*10^8)
 (7)
 > G:=[ piecewise( xS0, x/(1-DD)-E*Y               ),      piecewise( xS0, x/(1-DD)-E*Y               ),      piecewise( xS0, x/(1-DD)-E*Y               )    ]: plot(G, x=0..10^9, discont=true); fsolve~(G); SolveEquations~(G);
 (8)
 >

## Like mmcdara...

I'm not really sure what you are trying to achieve.

The attached will

1. Plot five more or less random graphs
2. Extract the data from these plots and build them into a matrix where the first column contains x-values: subsequent five columns contain y-values for each of the five graphs in (1) above
3. Export this matrix (in csv format - many other formats are possible)
 > restart;   with(plots):   with(plottools): # # Define five more or less random functions to plot #   v:= Vector       ( 5,         [1.1*t, 1.2*t, 1.3*t, 1.4*t, 1.5*t]       ): # # Plot the five functions above #   g11:= plot         ( [ seq             ( v[j],               j=1..5             )           ],           t = 1 .. 2.9,           color = [red, green, blue, black, purple],           numpoints = 200         ); # # Extract the plot data from the above in the form of # a matrix. # # First column is x-values: subsequent columns are the # y-values associated with each of the individual plots #   M1:= Matrix        ( [ getdata(g11)[1][3][..,1],            seq            ( getdata(g11)[j][3][..,2],              j=1..5            )          ],          scan=rows        ); # # Export this matrix to a csv file on my desktop. # # 1. Obviously the OP will have to change the #    location of the output file to something #    appropriate for his/her machine/OS # 2. Output will be in plain vanilla 'csv' format. #    Many other output formats are possible. #   ExportMatrix   ( "C:/Users/TomLeslie/Desktop/plotData.csv",     M1,     target=csv   );
 (1)
 >
 >

## Maybe a version issue...

The following works just fine in Maple 2019.1

```  restart;
interface(version);
with(Maplets[Elements]):
maplet1 := Maplet(["Hello World!", Button("OK", Shutdown())]):
Maplets[Display](maplet1);
```

## My \$0.02...

on why the PDE system in the worksheet Test_1.mw has got problems

1. Wrong number of BCs/ICs - you need five but you supply eight! This is simple to fix.
2. One of the PDEs in the system is diff(svt1(x, t), t) = 0. This means that svt1() cannot be a function of 't'. This can be fixed by removing this pde from the system, deleting the associated initial condition, and substituing svt1(x) for svt1(x,t) in the remaining equations
3. However one of the other PDEs in the system is diff(svt1(x, t), x) = 6342.941220/Lt1(x, t). Item (2) above shows that the left-hand-side of this equation is a function of 'x' only. Hence the right-hand-side must be a function of 'x' only. Hence Lt1(x,t) cannot be a function of 't'
4. Is this really what you expect????

See detailed diagnostics in the attached

 > restart:Digits:=30: Ts:=60*10^6*365*24*60*60: Tf:=2*Ts: phi0:=0.72: rho0:=646.8: E0:=1*10^9: nu0:=0.33: pc0:=4.0*10^6: mp:=1.0: pv0:=5.0*10^6: mvp:=0.8: g:=9.80665: Mp:=(6000*rho0)/Ts: lambda0:=E0*nu0/((1+nu0)*(1-2*nu0)): mu0:=E0/(2*(1+nu0)): K0:=lambda0+2/3*mu0: mus:=(1/16)*(-6*sqrt((K0-2*mu0)^2*phi0^2-(3*(K0-(8/9)*mu0))*(K0+2*mu0)*phi0+(9/4)*(K0+(8/9)*mu0)^2)-(6*(K0+2*mu0))*phi0+9*(K0-(8/9)*mu0))/(phi0-1): ks:=K0/(1-(1+3*K0/(4*mus))*phi0): eta:=1*10^9*10^6*365*24*60*60: a:=1.1545: A:=evalf(-2/3+a*sqrt(3)/3): B:=evalf(-1/3-a*sqrt(3)/3): Md:=t->Mp*t: phi:=1-(1-phi0)/Jir(x): pc:=pc0*(ln(phi)/ln(phi0))^mp: hp:=-(mp*(1-phi0)*pc)/(Jir(x)*phi*ln(phi)): pv:=pv0*(ln(phi)/ln(phi0))^mvp: hvp:=-(mvp*(1-phi0)*pv)/(Jir(x)*phi*ln(phi)): K:=4*ks*mus*(1-phi)/(3*ks*phi+4*mus): Kb:=-(1/3)*(3*ks+4*mus)/(Jir(x)*(3*ks+4*mus)-3*ks*(1-phi0)): mu:=mus*(1-phi)*(9*ks+8*mus)/(ks*(9+6*phi)+mus*(8+12*phi)): mub:=(-5/3)*(3*ks+4*mus)/(3*ks*(5*Jir(x)-2*(1-phi0))+4*mus*(5*Jir(x)-3*(1-phi0))): Gp1:=-(K+2*sqrt(3)/3*mu*a)/(K+a^2*mu+hp): Gp2:=-((sv(x)+2*sh(x))*Kb+a*sqrt(3)*(sv(x)-sh(x))*mub)/(K+a^2*mu+hp): Gvp1:=-(K+2*sqrt(3)/3*mu*a)/(K+a^2*mu+hp): Gvp2:=-((sv(x)+2*sh(x))*Kb+a*sqrt(3)*(sv(x)-sh(x))*mub)/(K+a^2*mu+hp): Gv1:=-(K+2*sqrt(3)/3*mu*a)/(K+a^2*mu+hvp): Gv2:=-((sv(x)+2*sh(x))*Kb+a*sqrt(3)*(sv(x)-sh(x))*mub)/(K+a^2*mu+hvp): Fhe:=K-2/3*mu: Fve:=K+4/3*mu: Fh:=((sv(x)+2*sh(x))*Kb-(sv(x)-sh(x))*mub): Fv:=((sv(x)+2*sh(x))*Kb+2*(sv(x)-sh(x))*mub): Fhp1:=Fhe+(K-sqrt(3)/3*mu*a)*Gp1: Fhp2:=(K-sqrt(3)/3*mu*a)*Gp2+Fh: Fvp1:=Fve+(K+2*sqrt(3)/3*mu*a)*Gp1: Fvp2:=(K+2*sqrt(3)/3*mu*a)*Gp2+Fv: Fhvp1:=Fhe+(K-sqrt(3)/3*mu*a)*Gvp1: Fhvp2:=(K-sqrt(3)/3*mu*a)*Gvp2+Fh: Fvvp1:=Fve+(K+2*sqrt(3)/3*mu*a)*Gvp1: Fvvp2:=(K+2*sqrt(3)/3*mu*a)*Gvp2+Fv: Fhv1:=Fhe+(K-sqrt(3)/3*mu*a)*Gv1: Fhv2:=(K-sqrt(3)/3*mu*a)*Gv2+Fh: Fvv1:=Fve+(K+2*sqrt(3)/3*mu*a)*Gv1: Fvv2:=(K+2*sqrt(3)/3*mu*a)*Gv2+Fv: He:=629.81: Hp:=1412.21: Hvp:=2403.25: Hs:=2798.10: Le:=0.997303781223145045267439740594: Lp:=0.446364534701243987: Lvp:=0.305960439144176588: Jire:=1: Jirp:=0.447352635989884784: Jirvp:=0.306391087534055584: she:=-1.97026501562798356531403389299*10^6: shp:=-7.27091982399678383000000*10^6: shvp:=-2.42633766270405289000000*10^7: sve:=-4.00023503172954239018303851002*10^6: svp:=-1.19729258539294481000000*10^7: svvp:=-2.98552061557406603000000*10^7: uh:=-2.25710511476647640150403675387*10^(-12): ue:=-8.54965365567908020218245625954*10^(-15): up:=-1.74701689029014799084132788157*10^(-12): uvp:=-4.45218466378321332440604021849*10^(-13):
 > # x to [Hs-He..Hs] unassign(sv,sh,L,Jir,u): sys_ode:= diff(sv(x),x)=rho0*g/L(x), diff(sh(x),x)=(K0-2/3*mu0)*diff(L(x),x)/L(x), diff(L(x),x)=rho0*g/(K0+4/3*mu0), diff(Jir(x),x)=0, diff(u(x),x)=-Mp*g/(K0+4/3*mu0): dvars:=indets({sys_ode},specfunc(anything,diff)): SYS:=solve({sys_ode},dvars): ics:=sv(Hs)=0,sh(Hs)=0,L(Hs)=1,Jir(Hs)=1,u(Hs)=uh: sol:=dsolve(SYS union {ics},numeric,output=listprocedure): sv1,sh1,L1,Jir1,u1:=op(subs(sol,[sv(x),sh(x),L(x),Jir(x),u(x)])): phi1:=x->1-(1-phi0)/Jir1(x): pc1:=x->pc0*(ln(phi1(x))/ln(phi0))^mp: pv1:=x->pv0*(ln(phi1(x))/ln(phi0))^mvp: rho1:=x->rho0/L1(x):
 > # x to [Hs-Hp..Hs-He] unassign(sv,sh,L,Jir,u): sys_ode:= diff(sv(x),x)=rho0*g/L(x), diff(sh(x),x)=Fhp1*diff(L(x),x)/L(x)+Fhp2*diff(Jir(x),x), Fvp1*diff(L(x),x)+Fvp2*L(x)*diff(Jir(x),x)=rho0*g, diff(Jir(x),x)/Jir(x)=-Gp1*diff(L(x),x)/L(x)-Gp2*diff(Jir(x),x), diff(u(x),x)=-Mp*g/(Fvp1+Fvp2*L(x)*diff(Jir(x),x)/diff(L(x),x)): dvars:=indets({sys_ode},specfunc(anything,diff)): SYS:=solve({sys_ode},dvars): ics:=sv(Hs-He)=sve,sh(Hs-He)=she,L(Hs-He)=Le,Jir(Hs-He)=Jire,u(Hs-He)=uh-ue: sol:=dsolve(SYS union {ics},numeric,output=listprocedure): sv2,sh2,L2,Jir2,u2:=op(subs(sol,[sv(x),sh(x),L(x),Jir(x),u(x)])): phi2:=x->1-(1-phi0)/Jir2(x): pc2:=x->pc0*(ln(phi2(x))/ln(phi0))^mp: pv2:=x->pv0*(ln(phi2(x))/ln(phi0))^mvp: rho2:=x->rho0/L2(x):
 > # x to [Hs-Hvp..Hs-Hp] unassign(sv,sh,L,Jir,u): sys_ode:= diff(sv(x),x)=rho0*g/L(x), diff(sh(x),x)=Fhvp1*diff(L(x),x)/L(x)+Fhvp2*diff(Jir(x),x), Fvvp1*diff(L(x),x)+Fvvp2*L(x)*diff(Jir(x),x)=rho0*g, diff(Jir(x),x)/Jir(x)=-Gvp1*diff(L(x),x)/L(x)-Gvp2*diff(Jir(x),x), diff(u(x),x)=-Mp*g/(Fvvp1+Fvvp2*L(x)*diff(Jir(x),x)/diff(L(x),x)): dvars:=indets({sys_ode},specfunc(anything,diff)): SYS:=solve({sys_ode},dvars): ics:=sv(Hs-Hp)=svp,sh(Hs-Hp)=shp,L(Hs-Hp)=Lp,Jir(Hs-Hp)=Jirp,u(Hs-Hp)=uh-ue-up: sol:=dsolve(SYS union {ics},numeric,output=listprocedure): sv3,sh3,L3,Jir3,u3:=op(subs(sol,[sv(x),sh(x),L(x),Jir(x),u(x)])); phi3:=x->1-(1-phi0)/Jir3(x): pc3:=x->pc0*(ln(phi3(x))/ln(phi0))^mp: pv3:=x->pv0*(ln(phi3(x))/ln(phi0))^mvp: rho3:=x->rho0/L3(x):
 (1)
 > # x to [0..Hs-Hvp] unassign(sv,sh,L,Jir,u): sys_ode:= diff(sv(x),x)=rho0*g/L(x), diff(sh(x),x)=Fhv1*diff(L(x),x)/L(x)+Fhv2*diff(Jir(x),x), Fvv1*diff(L(x),x)+Fvv2*L(x)*diff(Jir(x),x)=rho0*g, diff(Jir(x),x)/Jir(x)=-Gv1*diff(L(x),x)/L(x)-Gv2*diff(Jir(x),x), diff(u(x),x)=-Mp*g/(Fvv1+Fvv2*L(x)*diff(Jir(x),x)/diff(L(x),x)): dvars:=indets({sys_ode},specfunc(anything,diff)): SYS:=solve({sys_ode},dvars): ics:=sv(Hs-Hvp)=svvp,sh(Hs-Hvp)=shvp,L(Hs-Hvp)=Lvp,Jir(Hs-Hvp)=Jirvp,u(Hs-Hvp)=uh-ue-up-uvp: sol:=dsolve(SYS union {ics},numeric,output=listprocedure): sv4,sh4,L4,Jir4,u4:=op(subs(sol,[sv(x),sh(x),L(x),Jir(x),u(x)])): phi4:=x->1-(1-phi0)/Jir4(x): pc4:=x->pc0*(ln(phi4(x))/ln(phi0))^mp: pv4:=x->pv0*(ln(phi4(x))/ln(phi0))^mvp: rho4:=x->rho0/L4(x):
 > unassign('x','t'): Xp:=Hs-Hp: Xvp:=Hs-Hvp: phit1:=(x,t)->1-(1-phi0)/Jirt1(x,t): pvt1:=(x,t)->pv0*(ln(phit1(x,t))/ln(phi0))^mvp: fvpt1:=(x,t)->A*sht1(x,t)+B*svt1(x,t)-pvt1(x,t): Kt1:=(x,t)->4*ks*mus*(1-phit1(x,t))/(3*ks*phit1(x,t)+4*mus): Kbt1:=(x,t)->-(1/3)*(3*ks+4*mus)/(Jirt1(x,t)*(3*ks+4*mus)-3*ks*(1-phi0)): mut1:=(x,t)->mus*(1-phit1(x,t))*(9*ks+8*mus)/(ks*(9+6*phit1(x,t))+mus*(8+12*phit1(x,t))): mubt1:=(x,t)->(-5/3)*(3*ks+4*mus)/(3*ks*(5*Jirt1(x,t)-2*(1-phi0))+4*mus*(5*Jirt1(x,t)-3*(1-phi0))): Fhet1:=(x,t)->Kt1(x,t)-2/3*mut1(x,t): Fvet1:=(x,t)->Kt1(x,t)+4/3*mut1(x,t): Fht1:=(x,t)->((svt1(x,t)+2*sht1(x,t))*Kbt1(x,t)-(svt1(x,t)-sht1(x,t))*mubt1(x,t)): Fvt1:=(x,t)->((svt1(x,t)+2*sht1(x,t))*Kbt1(x,t)+2*(svt1(x,t)-sht1(x,t))*mubt1(x,t)): pde_sys_t1:={diff(Jirt1(x,t),[t])/Jirt1(x,t)=-fvpt1(x,t)/eta, diff(svt1(x,t),[t])=0, diff(svt1(x,t),[x])=rho0*g/Lt1(x,t), diff(sht1(x,t),[t])=Fhet1(x,t)*diff(Lt1(x,t),[t])/Lt1(x,t)-(Kt1(x,t)-sqrt(3)/3*mut1(x,t)*a)*diff(Jirt1(x,t),[t])/Jirt1(x,t)+Fht1(x,t)*diff(Lt1(x,t),[t])}: ics_t1:={Lt1(x,Ts)=L3(x),Lt1(Xp,t)=Lp,Jirt1(x,Ts)=Jir3(x),Jirt1(Xp,t)=Jirp,sht1(x,Ts)=sh3(x),sht1(Xp,t)=shp,svt1(x,Ts)=sv3(x),svt1(Xp,t)=svp}:
 > # # Check OP's code #   pds_t1:=pdsolve(pde_sys_t1,ics_t1,numeric,time=t,range=Xvp..Xp,output=listprocedure);
 > ############################################## #                                            # #     Diagnostic process                     # #                                            # ############################################## # # Eyeball how many derivatives of # what kind exist in the PDE system #   indets(pde_sys_t1, specfunc(diff));
 (2)
 > # # Output from the above command indicates that # the system needs *exactly* FIVE BCs/ICs. # # How many BCS/ICS is OP providing #   numelems(ics_t1);
 (3)
 > # # Need five BCs/ICs, and OP is providing eight. # So this isn't going to work. Might be able # to select five appropriate conditions from the # OP's eight??? # # The following will select the correct number of # of BCs/ICs of the correct type. From the output # of the indets() command above, the PDE system # needs the following conditions # #  1. Jirt1(x,t): IC (ie fixed t) #  2. Lt1(x,t):   IC (ie fixed t) #  3. sht1(x,t):  IC (ie fixed t) #  4. svt1(x,t):  IC (ie fixed t) #  5. svt1(x,t):  BC (ie fixed x) #   myIcs:=[ ics_t1[2], ics_t1[4], ics_t1[6], ics_t1[7], ics_t1[8]];
 (4)
 > # # Try for a solution with this reduced set of ICs/BCs #   pds_t1:=pdsolve(pde_sys_t1,myIcs,numeric,time=t,range=Xvp..Xp,output=listprocedure);
 > # # Never seen that error message before! - now what's wrong? # Examine the PDEs individually, for anything "anomalous". # # The third entry in the PDE system is interesting #    pde_sys_t1[3];
 (5)
 > # # If the derivative of svt1(x,t) with respect to 't' # is zero, then svt1(x,t) DOES NOT DEPEND on 't', which # can be illustrated by looking at the solution for # this equation on its own - which shows that it must be # a function only of 'x' #   pdsolve( pde_sys_t1[3]);
 (6)
 > # # One could(?) delete this PDE from the system, and maybe # substitute svt1(x) for svt1(x,t) in the remaining equations, # as well as removing the associated initial condition. # # But check out the fourth entry in the PDE system #    pde_sys_t1[4];
 (7)
 > # # Since svt1(x,t) is a function ONLY of 'x', the left-hand-side # of the above PDE is a function ONLY of 'x', and hence the # right-hand-side of this PDE is also a function ONLY of 'x', # and therefore Lt1(x,t) is a function ONLY of 'x'. # # This can be illustrated by producing general solutions for # the pair of equations formed by entries 3 and 4 in the PDE # system # # As I thought, both Lt1(x,t) and svt1(x, t) are functions # only of 'x' #   pdsolve( pde_sys_t1[3..4]);
 (8)
 > # # Again one ?could? substitute Lt1(x) for Lt1(x,t) everywhere # in the PDE system, but at least one of the remaining PDEs # depends on the derivative of Lt1() with respect to 't', so # quite a few terms in the remaining PDE system are now going # to be identically zero # # Probably not worth it - more likely that OP has some logical/ # conceptual problem with this PDE system
 >

## If...

you are interested in "playing" with random number generators, it is relatively easy to write your own package which will replicate the bahaviour of the in-built generators, at least for simple cases. In the attached I have provided a very simple package for dice-rolling (ie producing the numbers 1..6) which uses the same algorithm as the RandomTools[LinearCongruence] sub-package.

The worksheet shows that this simple package replicates the behaviour of the RandomTools[LinearCongruence] for the dice-rolling problem

NB Although it "works" - it is pretty basic, but it is only about 15 lines long! You could add *lots* of "bells and whistles" if you want

 > restart;   with(RandomTools[LinearCongruence]):   M1:=NewGenerator(range=1..6): # # Produce 100 numbers from this generator #   randInt1:=[seq( M1(), j=1..100)]; # # Get the state of the generator after the # above numbers have been produced #   GetState();
 (1)
 > # # Re-initialise the above random number generator #   SetState(state=1): # # Define another generator #   N1:=NewGenerator(range=1..6): # # Produce 100 numbers by alternating between the # two generators - just to see how "independent" # they actually are #   randInt2:=[seq( [M1(),N1()][], j=1..50)]; # # Verify that producing 100 numbers from a single # generator will give the same sequence as # alternating between two generators. Check the # state after the numbers have been produced: thi # should also be the same as before #   evalb(randInt1=randInt2);   GetState();
 (2)
 > # # Unload the above in-built PRNG #   unwith(RandomTools[LinearCongruence]):
 > # # Write code for a really simple package which replicates # the behaviour of the above "simple-case" PRNG #   myRand:= module()                  local mySeed, myGen:                  export NewGen, getSeed, setSeed:                  option package:                  mySeed:=1:                  NewGen:= proc()                                return(myGen):                           end proc:                  myGen:=proc()                              mySeed:=irem(427419669081*mySeed, 999999999989);                              return irem(mySeed, 6) + 1;                        end proc;                  getSeed:=proc()                                return mySeed;                           end proc;                  setSeed:=proc(n)                               mySeed:=n:                           end proc;            end module:   with(myRand):
 > # # Produce the first 100 numbers from this # "homebrew" PRNG #   M2:=NewGen():   randInt3:=[seq( M2(), j=1..100)]; # # Check that this sequences of random numbers is # exactly the same as that produced by the # "internal" generator above. Check the state # of the "homebrew" generator - should be the same # as that above #    evalb(randInt3=randInt1);    getSeed();
 (3)
 > # # Re-initialise the homebrew random number generator #   setSeed(1): # # Define another generator #   N2:=NewGen(): # # Produce 100 numbers by alternating between the # two generators - just to see how "independent" # they actually are #   randInt4:=[seq( [M2(),N2()][], j=1..50)]; # # Verify that produing 100 numbers from a single # generator will give the same sequence as # alternating between two generators. Check the # state after the numbers have been produced #   evalb(randInt3=randInt4);   getSeed();
 (4)
 >

## On reflection...

would't it be easier to use the geometry pacjkage for this?

As in the attached

 > with(geometry): point(A, [-5, 6]): point(B, [7, -3]): MakeSquare(sqr, [A, B, diagonal]): draw(sqr, axes=none);
 > detail(sqr);;
 (1)
 >

## Not sure why you have an error...

because the following works.

 > restart; interface(version); car_2som_opp := proc (U::list, V::list)  #construction d'un carré connaissant 2 sommets opposés                      local dist, eqCerU, eqCerV, r, sol, X, Y;                      dist := proc (M, N)                                   sqrt(add((M[i]-N[i])^2, i = 1 .. 2))                              end proc;                      r := dist(U, V)/sqrt(2);                      eqCerU := (x-U[1])^2+(y-U[2])^2 = r^2;                      eqCerV := (x-V[1])^2+(y-V[2])^2 = r^2;                      sol := solve([eqCerU, eqCerV], [x, y],allsolutions,explicit);                        map(allvalues,sol):                      X := [subs(op(sol[1]), x), subs(op(sol[1]), y)];                      Y := [subs(op(sol[2]), x), subs(op(sol[2]), y)];                      plots:-display(plot([U, X, V, Y, U],scaling = constrained, axes = none))                  end proc: car_2som_opp([-5,6],[7,-3]);
 >

## Not completely clear what you want...

and you will need at keast one boundary condition (which I have "guessed"), but maybe something in the attached will cover your requirement

 > restart; # # Define ODE using piecewise() #   ode:= diff(B[1](t),t)= piecewise                          ( t<=1000,                            kaC*(R-B[1](t))-kd1*B[1](t),                            -kd1*B[1](t)                          ); # # Add a "guesswork" boundary condition and # solve the ODE #   sol:=dsolve([ode, B[1](0)=1]); # # Plot the ODE solutions for values of # kd1 fromm 7e-03 to 7e-02 #   plot( [ seq           ( eval             ( rhs(sol),               [kaC=6e-02, kd1=j, R=1]             ),             j=7e-03..7e-02, 1e-03           )         ],         t=0..2000       ); # # Plot the ODE solutions for values of # kaC from 6e-02 to 6e-01 #   plot( [ seq           ( eval             ( rhs(sol),               [kaC=j, kd1=7e-03, R=1]             ),             j=6e-02..6e-01, 1e-02           )         ],         t=0..2000       ); # # Plot B[1](t) versus t and kd1, with kaC # as a parameter (produces mulitple surfaces) #   colors:=[red, green, blue, brown, black, yellow]:   plots:-display          ( [ seq              ( plot3d                ( eval                  ( rhs(sol),                    [kaC=6e-02+(j-1)*2e-02, R=1]                  ),                  t=0..2000,                  kd1=7e-03..7e-02,                  color=colors[j]                ),                j=1..6              )            ]          );
 >

## The command...

numbers := [random[empirical[0.5, 0.5]](N)]

uses the 'stats' package, which has been 'deprecated' since Maple version 9.5 in 2004. You cannot access these commands using 'with(Statistics)'. You have to use 'with(stats)' , as shown in the attached - which actually "works".

Note that you really shouldn't be using such deprecated packages without a very very good reason. Since all this command actually does is produce a list containing the integers 1,2 at random - there really doesn't seem to be any sensible reason for continuing to use it!

Also

1. in the final for loop: when you are just adding a finite list of "known" numbers, it is much more efficient to use add() rather than sum()
2. the procedure Wpath() is defined but never used - deliberate?

 >
 >
 >
 >
 >
 >
 >
 >
 (1)
 >

## What exactly is wrong?...

Somehow, the "code" in your file Matrix.mw appears to be "Maple Output" which is odd

If I convert everything to 1-D Maple Input, then everything seems to execute correctly - see the attached

 > restart;   interface(rtablesize = 10);
 > with(LinearAlgebra);
 > A := 8:   B := 5:   q := 0.4:   p := 0.2:   r := 1 - p - q:   dimP := A + B + 1:
 > P := Matrix(dimP, dimP, [0 \$ dimP*dimP]):   P[1, 1] := 1:   P[1, 2] := 0:   P[dimP, dimP] := 1:   P[dimP, dimP - 1] := 0:   for i from 2 to dimP - 1 do       P[i, i - 1] := q;       P[i, i] := r;       P[i, i + 1] := p;   end do:
 > p0 := Matrix(dimP, 1, [0 \$ dimP]):   p0[A + 1, 1] := 1:   pV[0] := p0:   PT := Transpose(P):
 > for n to 200 do       pV[n] := PT . (pV[n - 1]):   end do:
 > map(x -> evalf(x, 3), Transpose(pV[5]));
 (1)
 >

## Difficult to tell...

what you are getting wrong unless you upload the relevant worksheet using the big green up-arrow in the toolbar.

The attached will produce the plot you want

 > plot(x^2-9, x=-10..10);
 >

## One way...

is to get rid of the integral term, by differentiating the pde with respect to 't'.

The drawback with this approach is hte at it increases the degree of the PDE (with respect to 't') so an additional initial condition is needed. I tried a few possibilities for D[2](u)(x,0). These didn't make a huge difference to the overall "form" of the solution, whihc always gets *seriously* huge for t>~1.
You might want to experiment with the following, and in particular with the "extra" initial condition D[2](u)(x,0=1, which I inserted. (Changing the rhs of this initial condition doesn't make a huge difference to the overall form!)

 > interface(rtablesize=10): f:=(x^3+t^2*x^2-t^2*x+4*t*x-2*t+1)*exp(x*t)-x+1; pde:=diff(u(x,t),t)+diff(u(x,t),x,x)+u(x,t)+int(u(x,tau),tau=0..t)=f; pde:=diff(pde,t); bounds:= u(0,t)=0, u(1,t)=0, u(x,0)=x*(x-1), D[2](u)(x,0)=1; psol:=pdsolve(pde, [bounds], numeric); psol:-plot3d( u(x,t), x=0..1, t=0..0.1);
 >

## Something like...

the attached -maybe?

 > restart;   interface(rtablesize=12);   M2:=Matrix(2, 12, (i,j)->`if`( i=1,                                  `if`(type(j, odd), U[(j+1)/2],0),                                  `if`(type(j, even), U[j/2], 0)                                )             );   Matrix( [ [seq( D[1](M2[1,j])(x,y) + 0*M2[2,j],          j =1..12)],             [seq( 0*(M2[1,j])(x,y)   + D[2](M2[2,j])(x,y), j =1..12)],             [seq( D[2](M2[1,j])(x,y) + D[1](M2[2,j])(x,y), j =1..12)]           ]         ):   convert~([%], diff);
 (1)
 >