## 1585 Reputation

19 years, 261 days

## Use a procedure....

There was an error in your code, I hope you can spot it.

Use this:

```F1 := proc(Mh := 1)
local m, H, b, a, V, k, p, j, v, chi, appr, u_appr, u_appr_1;
m := 10;
H := 1;
b := 0.02;
a := 0.05;
V := array(0 .. m);
V[0] := 1 - exp(-t);
for k to m do if k = 1 then chi := 0 else chi := 1 end if;
p := 0;
for j from 0 to k - 1 do p := p + V[k - 1 - j]*(diff(V[j], `\$`(t, 2)))
- (diff(V[k - 1 - j], t))*(diff(V[j], t)) - a*(2*(diff(V[k - 1 - j], t)
)*(diff(V[j], `\$`(t, 3)))
- (diff(V[k - 1 - j], `\$`(t, 2)))*(diff(V[j], `\$`(t, 2)))
- V[k - 1 - j]*(diff(V[j], `\$`(t, 4))))
end do;
p := (p + (diff(V[k - 1], `\$`(t, 3)))
- b*((diff(V[k - 1], `\$`(t, 2))) + t*(diff(V[k - 1], `\$`(t, 3))))
- Mh*(diff(V[k - 1], t)))*h*H;
p := factor(p);
V[k] := -(int(p, t)) + 0.5*exp(t)*(int(exp(-t)*p, t))
+ 0.5*exp(-t)*(int(exp(t)*p, t)) + chi*V[k - 1] + C1 + C3*exp(-t);
v := unapply(V[k], t);
V[k] := frontend(expand, [V[k]]);
V[k] := subs(C3 = solve(eval(subs(t = 0, diff(V[k], t))), C3), V[k]);
V[k] := frontend(expand, [V[k]]);
V[k] := subs(C1 = solve(eval(subs(t = 0, -V[k] - (diff(V[k], t)))), C1),
V[k]);
end do;
appr := 0;
for k from 0 to m do appr := appr + V[k] end do;
u_appr := unapply(appr, h, t);
u_appr_1 := unapply(diff(u_appr(h, t), t), h, t);
evalf(u_appr_1(-0.4, t));
plot([u_appr_1(-0.4, t)], t = 0 .. 4, 0 .. 1.2, color = [black], axes = frame
);
end proc;

plots[display](F1(2), F1(3));
```

Thumb if You like.

## Just:...

eval(ode, y=exp);

Thumb if you like.

## Correct the wrong definitions....

Like this:

```bix:= [seq(1..4,0.1)];  NN := nops(bix);
prx:=[seq(1..2,0.1)];  NN1 := nops(prx);
for i  from 1 to NN do
for j from 1 to NN1 do
R := dsolve(eval({bc, eq1,eq2}, {bi=bix[i],pr=prx[j]}), fcns, type = numeric, method = bvp[midrich], maxmesh=2400):
X1[i,j]:=rhs(-R(0)[3]):
end do:
end do:
```

Thumb if You like.

## This question is so frequent....

See my answer to an older query:

```function_coeffs := proc(A, v::set(name))
local S, T;
S := indets(A, {function});
S := select(has, S, v);
T := {Non(map(identical, S))};
frontend(proc(A, S) local V; [coeffs](collect(A, S, distributed), S, 'V'), [V] end proc, [A, S union v], [T, {}])
end proc:

eq := 1.1*cos(t) + 0.8*sin(t) + t^2*x + x;

function_coeffs(eq, {t});
```

Thumb if You like.

## Just some formatting......

Use this:

```sol := solve([seq(E[i] = 0, i = 1 .. 5)], [S, H, r, E, P]):
sol := collect(evala(sol), [RootOf, mu, a, b], factor):
alias(Q = indets(sol, RootOf)[1]):
subs(Q = q, sol);
collect(collect((evala@Norm)(q - Q), q, normal), [q, nu, a, b] ,factor);
```

Also notice that you are using E as a table (defining E[1..5]) and as variable in equations. This will fail you eventually.

Thumb if You like.

## Primpart....

Lookup primpart, content. These can be used to separate variables from coefficients of a monomial. Likewise you could use coeff.

Example

```K := proc(monomials :: set, V :: {set, list})
description "Test if the set is a set of multiples of the same monomial in variables V";
map((x -> primpart(collect(x, V, distributed), V)), monomials);
if not type(`%`, set(monomial)) then error "Input is not a set of monomials"; fi;
evalb(nops(`%`) = 1);
end;

K({3*(a+b)*(a-b)+3*b^2, q*a^2}, {a, b});
K({3*(a+b)*(a-b)+3*b^2, q*a^2}, {q});
```

Thumb if You like.

## Like this:...

```S := remove(type, remove(has, indets(B, function), diff), trig);
frontend(proc(E, S) local V; [coeffs(collect(E, S, distributed), S, 'V')], [V]; end, [B, S], [{Non(function)}, {}]):
V := collect([%], [diff], distributed, factor);
```

or like this:

```S := remove(type, indets(B, function), trig);
frontend(proc(E, S) local V; [coeffs(collect(E, S, distributed), S, 'V')], [V]; end, [B, S], [{Non(function)}, {}]):
V := collect([%], [diff], distributed, factor);
W := ListTools[Classify](proc(E) op(0, convert(E[2], D)); indets(%, name); end, zip(() -> [args], op(V))):
W := map(proc(S) [op](S); map2(map2, op, [1,2], %) end, W);
```

Thumb if You like.

## Lookup ?plots[pointplot3d]...

## Using reals....

Comonly you would use evalc under assuming real.

For example:
evalc(%) assuming real, R+r*cos(theta)>0, r>0;

Thumb if You like.

## Roots....

The command roots was specifically designed for finding rational roots:

```P:=`*`('(x*rand(-99..99)()+rand(-99..99)())^op(rand(1..4)(),[1\$3,2])'\$8)*
`*`('(x^2*rand(-99..99)()+rand(-99..99)())^op(rand(1..4)(),[1\$3,2])'\$8);
expand(P);
roots(expand(P));
```

Thumb if You like.

## Classify....

Note: the use of has can be weak.

```seq(i = nops(ListTools[Classify](has, Set, i)[true]), i = (`union`@op)(Set));

seq(i = nops(ListTools[Classify]((s,i)->member(i,s), Set, i)[true]), i = (`union`@op)(Set));
```

This one looks cleanest:

```map((i->i) = ListTools[Occurrences], ((`union`@op) , [op])(Set), member);
```

Without ListTools:

```F:=(Set::{set,list}({set,list})) -> map((i->i) = (i-> (nops@select)(has, Set, i)), ({op}@map)(op, Set));
F(Set);

F:=(Set::{set,list}({set,list})) -> map((i->i) = (i-> (nops@select)((s,i)->member(i,s), Set, i)), ({op}@map)(op, Set));
F(Set);
```

Thumb if You like.

## Use simplify....

This is the way to these substitutions, although I do not approve of using eval:

```expression_ABC:=eval(expression,[x^(-(3*c)/2 + a/2 + 1/2)*(c + a + 1)*b^(-(3*c + a + 1)/(2*c))=A,c*x^(a/2 + 1/2 - c/2)*b^(-(c + a + 1)/(2*c))=B]);
simplify(%, [((a + 1)*(c + a + 1)*(2*c + a + 1))=C]);
```

Thumb if You like.

There was a typo in your code.

MapleQuestions_227394.mw

Thumb if You like.

## Modules....

I would expect LCState is a variable that belongs the module containing the GetState proc. That would be one of the modules:  RandomTools:-LinearCongruence and RandomTools.