Question:How I can convert my maple code to matlab?

Question:How I can convert my maple code to matlab?

Maple 2018

How I can convert my maple code to matlab?.

 > restart;

 > fdsolve := proc({alpha:=NULL, t0:=NULL,                  t1:=NULL, x0:=NULL, y0:=NULL,                  N:=NULL}, params)     local t, h, h1, h2, c, b, x, y, L, n, l, X, Y, f, g;     eval(F(t,x,y), params);     f := unapply(%, [t,x,y]);     eval(G(t,x,y), params);     g := unapply(%, [t,x,y]);     L := floor(1/alpha);     h := (t1 - t0)/N;     h1 := h^alpha/GAMMA(alpha+1);     h2 := h^alpha/GAMMA(alpha+2);     c := (i,n) ->         `if`(i=0,             (n-1)^(alpha+1) - (n-1-alpha)*n^alpha,             (n-i+1)^(alpha+1) + (n-i-1)^(alpha+1) - 2*(n-i)^(alpha+1));     b := (i,n) -> (n-i)^alpha - (n-1-i)^alpha;     t := Array(0..N, i-> (1-i/N)*t0 + i/N*t1, datatype=float[8]);     x[0], y[0] := x0, y0;     for n from 0 to N-1 do         X[0], Y[0] :=             x[0] + h1*add(b(i,n+1)*f(t[i],x[i],y[i]), i=0..n),             y[0] + h1*add(b(i,n+1)*g(t[i],x[i],y[i]), i=0..n);         for l from 1 to L do             X[l], Y[l] :=                 x[0] + h2*add(c(i,n+1)*f(t[i],x[i],y[i]), i=0..n)                      + h2*f(t[n+1], X[l-1], Y[l-1]),                 y[0] + h2*add(c(i,n+1)*g(t[i],x[i],y[i]), i=0..n)                      + h2*g(t[n+1], X[l-1], Y[l-1]);         end do;         x[n+1], y[n+1] := X[L], Y[L];         #printf("y[%d]=%a\n", n+1, y[n+1]);     end do;     return Array(0..N, i -> [t[i], x[i], y[i]]); end proc:
 >
 > F := (t,x,y) -> r*x*(1-x/k) - beta*x*y/(a+x^2); G := (t,x,y) -> mu*beta*x*y/(a+x^2) - d*y - eta*x*y;
 (1)
 > params := { r=0.05, a=0.8, mu=0.8, d=0.24,             eta=0.01, beta=0.6, k=1.6 };
 (2)
 > T := 300.0;  # time interval: 0 < t < T q := 100;    # divide the time interval into q subintervals
 (3)

We produce several solutions starting from various initial data .

This reproduces the phase diagran in the cited article's Figure 2  (alpha=0.98):

 > sol := fdsolve(alpha=0.98, t0=0.0, t1=T, x0=2.5, y0=0.14, N=q, params): p1 := plot([seq([sol[i][2], sol[i][3]], i=0..q)])
 >
 >