adams_solver.jl

function fun_linear_system(t, x, A)
  return A*x;
end

nx = 4;
eigens = diagm(-1.+(1:nx)/nx/4);
for i=1:2*floor(nx/2)
  eigens[1,i]=0.4^(i/2);
  eigens[i,1]=-0.4^(i/2);
end
T=randn(nx,nx);
A = T*eigens\T;
fun = (t,x)->fun_linear_system(t,x,A)

T=15
N=1000
α=0.9
y0 = rand(nx,1);

function solve_adams(fun, α, y0, T=30, N=3000)
  h = T/N;
  a = z-> (z+1)^(α+1) + (z-1)^(α+1) - 2*z^(α+1);
  b = z-> z^α - (z-1)^α;
  y = Array(Float64, length(y0), N)
  y[:,1]=y0;
  for n=1:N-1
    s = 0; for j=0:n-1 s += b(n+1-j) * fun(n*h, y[:,j+1]) end
    p = y0 + h^α*s/gamma(α)
    s=0; for j=1:n-1 s += a(n+1-j) * fun(n*h, y[:,j+1]) end
    y[:,n+1]  = y0 + (h^α/gamma(α+2))* (
            fun(n*h, p) + s  +
             ((n-1)^(α+1)- (n-1-α)*n^α)*fun(0,y[:,1])
            );
  end
  y
end

y = solve_adams(fun, α, y0, T, N)
plot(h*(1:N),y')