%%Ref: num-ht5/Data/2023-07-14-nogauge-hkn-sys.rr
xmax=200;
%xmax=300;
%xmax=400;
%xmax=500;
%xmax=1000;
A=chebop([1,xmax]);
yy=chebfun('y',[1,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))];

A.lbc= @(z0,z1,z2,z3) [z0-0.0781014; z1-0.0509628];
% y=200
A.rbc= @(z0,z1,z2,z3) [z0-2.253978898520279*10^9; z1-1.531999507742409*10^8];
% y=300  It seems to work.
%A.rbc= @(z0,z1,z2,z3) [z0-1.034983054985926*10^12; z1- 5.777933784872538*10^10];
% y=400  It makes a mild oscillation.
%A.rbc = @(z0,z1,z2,z3) [z0-1.860971664352828*10^14; z1-9.03084277029875*10^12];
% y=500  It does not work.
%A.rbc= @(z0,z1,z2,z3)  [z0-1.830333399187609*10^16; z1- 7.965451393048757*10^14];
%y=1000,  value at y=4 is 10^4! It does not work.
%A.rbc= @(z0,z1,z2,z3) [z0-1.322069561001868*10^24;   z1-4.096770343810229*10^22];

S=A\0;
%Compare, e.g., S{1}(4) and  myhkn(4)
%y2023_07_14_nogauge_hkn*.m のスクリプトファイルにして、コマンドwindow からファイル名で実行. 開いてから実行でもいい. 関数のみのファイルは起動時に自動で実行される.
%plot(S{1},'.-'); e.g., xlim([2,5]); xlim auto
disp(myhkn(4))
disp(S{1}(4))
xlim auto ; plot(S{1})
function hval = myhkn(y)
  xrange=1:1:40;
  yval=[0.0781014,0.140353,0.229823,0.354688,0.525001,0.753032,1.05367,1.44488,1.94825,2.58958,3.39957,4.41466,5.67789,7.23995,9.16035,11.5088,14.3666,17.8284,22.0043,27.0217,33.028,40.1932,48.7132,58.8128,70.7503,84.8212,101.363,120.763,143.458,169.951,200.808,236.675,278.282,326.456,382.132,446.364,520.342,605.407,703.065,815.011];
  hval=interp1(xrange,yval,y);
end