See this help page  rtable_indexing

From Fig. 2 it is clearly seen that  f(0)=1  and  D(f)(0)<0 , which contradicts your initial conditions.

 Jayaprakash J   For movable boundary conditions you can use the  Explore  command.

Download BVP_new.mw

restart;
  rnge:=-1..1:
  fcns:=[[[0,0]]$5, x$10, [x,x^2]$10, [x,x^2,(1-x)]$10]:
  plots:-display( [ seq
                    ( plot
                      ( fcns[j], x=rnge),
                      j=1..35
                    )
                  ],
                  insequence=true, scaling=constrained
                );

Here is a more complicated option:

restart;
with(plots):
A:=animate(plot, [x, x=-1..a, color=red, thickness=2], a=-1..1):
B:=animate(plot, [x^2, x=-1..a, thickness=2, color=green], a=-1..1):
C:=animate(plot, [1-x, x=-1..a, color=blue, thickness=2], a=-1..1):
display([A, display(op([1,-1,1],A),B), display(op([1,-1,1],A),op([1,-1,1],B),C)], insequence, scaling=constrained);

                

See this post   https://mapleprimes.com/posts/207840-Combinations-Of-Multiple-Animations

Edit.

Download orthogonal-new1.mw

 I think the easiest way to create the same animation is the using a procedure that creates an one frame, and then using the plots:-animate command:

restart;

OneFrame:=proc(a)
local f, Point, Points1, Points2, Curve;
uses plots;
f:=x->x^2-x+2;
Point:=plot([[2,f(2)]],style=point, color=black, symbol=solidcircle, symbolsize=30);
Points1:=plot([[2-a,f(2-a)],[2+a,f(2+a)]], x=-3..5, style=point, color=red, symbol=solidbox, symbolsize=30);
Points2:=plot([[2-a,0],[2+a,0],[0,f(2-a)],[0,f(2+a)]], x=-3..5, style=point, color=red, symbol=diamond, symbolsize=30);
Curve:=plot(f, -3..5, 0..8, color=blue, thickness=3);
display(Point,Points1,Points2,Curve, gridlines, scaling=constrained);
end proc:

plots:-animate(OneFrame, [a], a=1..0.05);

See  https://en.wikipedia.org/wiki/MapleSim

This can be done in many ways. Here is a simple procedure with a for-loop that solves the problem:

P:=proc(A::Matrix, j::posint, a::numeric)
local i, m, n;
m,n:=upperbound(A);
for i from 1 to m do
if A[i,j]>a then return i fi;
od;
end proc:

Example:
A:=<1,2,1.1; 3,4,2.1; 5,6,3.1>;
P(A,3,2);

To draw this circle centered at point  H  and radius  R , you do not need to know any angles. Simply write the appropriate parametric equations and use  plot3d  command:

restart; 
with(geom3d):
 a := 3:
 b := 4:
 h := 5:
 point(A, 0, 0, 0):
 point(B, a, 0, 0):
 point(C, a, b, 0):
 point(DD, 0, b, 0):
 point(S, 0, 0, h):
 sphere(s, [A, B, C, S], [x,y,z], 'centername' = m); 
detail(s); 
plane(p, [S, A, B], [x, y, z]); 
n:=coordinates(projection(H, m, p)); 
R := distance(H, S);
S1:=plot3d(n+r*~([cos(t),0,sin(t)]), r=0..R, t=0..2*Pi, style=surface,  color=pink):
S2:=plots:-spacecurve(n+[R*cos(t),0,R*sin(t)], t=0..2*Pi, color="VioletRed", thickness=2):
S3:=plottools:-sphere([3/2, 2, 5/2], 5*sqrt(2)*(1/2), color="LightGreen", transparency=0.7):
plots:-display(S1,S2,S3, axes=normal, orientation=[20,70], lightmodel=light1);

T := (f,k) -> 1/k!*diff(f,[x$k]);

Examples of use:

T(sin(3*x), 5);
T(sin(3*x), k);

In my opinion, Maple simplifies the expression quite well. But what is your main goal? If you want a simple way to calculate this expression for specific parameters values  a, b, c, d  and a specific  n , then here is a simple procedure for this:

restart;
Sol:=proc(a,b,c,d,n0)
simplify(rsolve({f(n+1)-((a-1)*(b-1)+(c+d))*f(n)-c*d*f(n-1)=0, f(0)=a,f(1)=a*(b - 1)+c},f));
simplify(eval(%,n=n0));
end proc:	

Examples of use:

Sol(2,1,4,3,100); ``;
seq(Sol(2,1,4,3,k), k=1..20); ``;
Sol(2,1,4,3,n);

f := x -> 3*x + 2:
g := D(f):
g(3);

P:=proc(p,q,r)
local pq, pr, v, cv, T;
uses LinearAlgebra, plots, plottools;
pq:=arrow(p,q-p, width=0.1, color=red);
pr:=arrow(p,r-p, width=0.1, color=red);
v:=CrossProduct(convert(q-p,Vector),convert(r-p,Vector));
cv:=arrow(p,v, width=0.1, color=green);
T:=polygon([p,q,r], color=yellow);
display(pq,pr,cv,T, scaling=constrained, axes=normal, lightmodel=light1);
end proc:

Example of use:

P([3,3,2], [3,4,5], [4,3,4]);

If one solution is enough, then you can do so

restart;
FourSquares:=proc(n)
local a, b, c, d;
for a from 1 to floor(sqrt(n/4)) do
for b from a to floor(sqrt((n-a^2)/3)) do
for c from b to floor(sqrt((n-a^2-b^2)/2)) do
for d from c to floor(sqrt(n-a^2-b^2-c^2)) do
if n=a^2+b^2+c^2+d^2 then return [a,b,c,d] fi;
od: od: od: od:
end proc:

Example of use:

FourSquares(100000);
                                               [4, 8, 8, 316]

It is not recommended to use some variable and indexed variable with the same name in the same code. Therefore, I replaced, for example, t[n]  by  t||n  and so on. Some other corrections and additions were also made. The resulting approximate solution is quite accurate.

restart;
f:=t->y(t)-t^2+1;
eqn:=diff(y(t),t)=f(t):
ex:=dsolve({eqn,y(0)=0.5},y(t));
t||0:=0: w||0:=0.5: h:=0.2: ex||0:=0.5: e||0:=0:
Dy||0:=eval(f(t),[t=t||0,y(t)=w||0]);

for n from 1 to 10 do
t||n:=n*h; ex||n:=eval(rhs(ex),t=t||n);
w||n:=w||(n-1)+h*Dy||(n-1)+h^2/2!*(Dy||(n-1)-2*t||(n-1));
e||n:=abs(ex||n-w||n);
Dy||n:=eval(f(t),[t=t||n,y(t)=w||n]);
od:

printf(" i | t[i] |(Taylor)w[i] |(exact)y[i] |Error | \n ");
for i from 0 to 10 do
printf("%2.2f| %5.2f  | %5.6f| %5.6f  |  %5.6f | \n", i, t||i, w||i ,ex||i,e||i) ;
od;

A:=plot(rhs(ex), t=0..2, color=red):
B:=plot([seq([h*k,w||k], k=0..10)], style=point,color=blue):
plots:-display(A,B, legend=["Exact solution","Approximate solution"], view=0..6);

The final results as a table and a plot: