%Ref: num-ht5/Data/2023-07-14-nogauge-hkn-sys.rr
%%Ref: num-ht5/Data/2023-07-25-hkn-nogauge-int.m (Mathematica)
xstart=10^4;
xmax=10^4+40;
mag=1;
magexp=82;
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))];

r=1+(0.004-rand(1,4)*0.008)
%r=1+(0.05-rand(1,4)/10);
%r=[1,1,1,1];
lz0 = r(1)*6.8650347328000341103*10^(82-magexp);
lz1 = r(2)*6.8164613890282423426*10^(80-magexp);
rz0 = r(3)*1.0208464183481783696*10^(83-magexp);
rz1 = r(4)*1.0116158445741086178*10^(81-magexp);
%y=10^4
A.lbc= @(z0,z1,z2,z3) [z0-mag*lz0; z1-mag*lz1];
%y=10^4+40. 
A.rbc = @(z0,z1,z2,z3) [z0-mag*rz0;  z1-mag*rz1];

S=A\0;
%plot(S{1},'.-'); e.g., xlim([2,5]); xlim auto
hold on
plot(S{1},'.-')
[lz0,rz0]