%%Ref: num-ht5/Data/2023-07-14-nogauge-hkn-sys.rr
%%Ref: num-ht5/Data/2023-07-25-hkn-nogauge-int.m (Mathematica)
%%Note: how far is it trustful?
xstart=10^8;
magexp=8678;

%modify the following lines.
% values of H^k_n and diff(H^k_n,y) at y=xmax
xmax = 100200000;
rz0 = 3.8569183252984497846*10^(8687-magexp);
rz1 = 3.8527786677379301901*10^(8683-magexp);
% loosing 3 digits accuracy at 10^8+200
%xmax = 100300000;
%rz0 = 8.3838241945947072011*10^(8691-magexp);
%rz1 = 8.3706501601553094774*10^(8687-magexp);

mag=1;
A=chebop([xstart,xmax]);
A.op= @(y,z0,z1,z2,z3) [diff(z0)-(z1); diff(z1)-(z2); diff(z2)-(z3); diff(z3)-(((y^2-4*y)*z3+(15*y-2)*z2+(-y+13)*z1-11*z0)/(y^2))];

%y=10^8
A.lbc= @(z0,z1,z2,z3) [z0-mag*8.04155070290508114583167116889111540523256835856; z1-mag*8.0409479557121269329004693340582168507703767*10^(8674-magexp)];
%y=xmax
A.rbc = @(z0,z1,z2,z3)  [z0-mag*rz0;z1-mag*rz1];

S=A\0;
%plot(log(S{1}),'.-'); e.g., xlim([2,5]); xlim auto
xlim auto
plot(S{1},'.-')
disp("S{1}(10^8+100) is 8.1223?")
S{1}(10^8+100)
S{1}(xmax-10^5)
