Numerical Solution To The 
Spreaded Given Problems For 
Linear Ordinary Differential Equations

Differential Equation Solver

(*:Source:
	1. A . Nersessian .  Generalization of the Euler-Maclaurin Formula ans Some
Applications. Theory of functions and applications. Collection of
 Works Dedicated to the Memory of Mkhitar M .  Djrbashian, "Louys" Publishing House, 
    1995, 133 - 138. 
	 *)

	  

	\( (*\(:\)\(Limitations : 
	
            Best\ if\ the\ coefficients\ of\ the\ equation\ right\ side\ and\ possible\ integrands\ are\ rather\ \(\(smooth\)\(.\)\)\)
	*) \)

BeginPackage["NumericalMath`NLinearDSolve`"];
Off[
        General::spell];
Off[General::spell1];
Off[
        Break::nofwd];
Off[NIntegrate::"inum"];
Off[
        NIntegrate::"nlim"];
Off[
        NIntegrate::"vars"];
\(NLinearDSolve::usage = "NLinearDSolve[eqns, y, {x, xmin, xmax}] finds a numerical solution to the initial or boundary        problems for ordinary linear differential equations for the function y with the independent variable x in        the range xmin to xmax.";\)
NLinearDSolve::ndode = "Input is not an \
ordinary linear  differential equation.";
\(NLinearDSolve::ndnco = "The number of conditions is not equal to the total differential        order of the equation.";\)
\(NLinearDSolve::fcabsun = "NLinearDSolve failed to converge to prescribed precision after using `1` points of the grid.        Current Uniform-norm (with weight->`4`) of the solution equals about `2`. Absolute Uniform-error (with weight->`4`) equals about `3`. For more        precision use more MaxGridPoints.";\)
NLinearDSolve::fcrelun = "NLinearDSolve failed to converge to prescribed \
precision after using `1` points of the grid.        Current Uniform-norm  \
(with weight->`4`) of the solution equals about `2`. Relative Uniform-error \
(with weight->`4`) equals about `3`. For more        precision use more \
MaxGridPoints.";
NLinearDSolve::fcabsL2 = "NLinearDSolve \
failed to converge to prescribed precision after using `1` points of the \
grid.        Current L2-norm (with weight->`4`) of the solution equals about \
`2`. Absolute L2-error (with weight->`4`) equals about `3`. For more \
precision use more MaxGridPoints.";
\(NLinearDSolve::fcrelL2 = "NLinearDSolve failed to converge to prescribed precision after using `1` points of the grid.        Current L2-norm (with weight->`4`) of the solution equals about `2`. Relative L2-error (with weight->`4`) equals about `3`. For more precision use        more MaxGridPoints.";\)
NLinearDSolve::unknnopt = "Unknown Option `1`.";
NLinearDSolve:: 
          errmes = "Unknown option for Error. Now NLinearDSolve works with \
Error->Relative.";
NLinearDSolve::symsol = "Unknown option for \
SymbolicSolution. Now NLinearDSolve works with \
SymbolicSolution->False.";
NLinearDSolve::precinf = "Unknown \
option for PrecisionInformation. Now NLinearDSolve works with \
PrecisionInformation->False.";
NLinearDSolve::nrmes = "Unknown \
option for SolutionNorm. Now NLinearDSolve works with \
SolutionNorm->L2.";
NLinearDSolve::comp = "Unknown option for \
Compiled. Now NLinearDSolve works with Compiled->True.";
NLinearDSolve::prg = "Unknown option for PrecisionGoal. Now NLinearDSolve \
works with PrecisionGoal->5.";
NLinearDSolve::mxgpmes = "Unknown option for MaxGridPoints. \
Now NLinearDSolve works with MaxGridPoints->200.";

Options[NLinearDSolve] = {SymbolicSolution -> False, 
            PrecisionInformation -> False, 
            PrecisionGoal -> Automatic, MaxGridPoints -> Automatic, 
            SolutionNorm -> L2, Error -> Relative, 
            SolutionNorm -> Uniform, Error -> Absolute, 
            Weight -> Automatic, InterpolationOrder -> 6, 
            InterpolationPoints -> 200, MaxRecursion -> 6, 
            WorkingPrecision -> $MachinePrecision, 
            AccuracyGoal -> $MachinePrecision - 10, 
            InterpolationPrecision -> $MachinePrecision - 10};

Begin["`Private`"];

RecognizeDif[eq_, y_, x_] := 
          Module[{ord, eq3, eq4, t1, t2, t3, t4, pp, t5, t6, rr, p, r, t7, 
              t8}, 
ord = 
              Max[Cases[eq, Derivative[_][_][_], 
                      Infinity] /. Derivative[n_][_][_] -> 
                      n] /. -Infinity -> 0; 
If[
              Positive[ord] && IntegerQ[ord], , 
              Message[NLinearDSolve::ndode]; 
              Break[]]; 
eq3 = 
              eq[[1]] - eq[[2]]; 
eq4 = 
              ReplaceRepeated[eq3, 
                k1_.  Derivative[i_][y][x] + 
                    k2_.  Derivative[i_][y][x] -> (k1 + 
                        k2) Derivative[i][y][
                      x]]; 
ClearAll[
              a]; 
t1 = 
              Cases[eq4 + 
                  a,   k_.  Derivative[i_][y][
                    x]]; 
t2 = 
              Reverse[Table[
                  If[Length[t1[[i]]] > 1, 
                    Drop[t1[[i]], {Length[t1[[i]]]}], 1], {i, 1, 
                    Length[t1]}]]; 
t3 = 
              Cases[eq4, k_.  y[x]]; 
t4 = 
              If[Length[t3] == 0, 0, 
                Apply[Plus, 
                  Table[If[Length[t3[[i]]] > 1, 
                      Drop[t3[[i]], {Length[t3[[i]]]}], 1], {i, 1, 
                      Length[t3]}]]]; 
pp = 
              Append[t2, t4]; 
t5 = 
              Position[eq4, y[x]]; 
t6 = 
              Position[
                eq4, Derivative[i_][y][
                  x]]; 
rr = -Take[
                  eq4, {1, (If[Length[t5] == 0, t6, t5] // First // 
                          First) - 1}]; 
p = 
              pp/pp[[1]]; 
r = 
              rr/pp[[1]]; 
t7 = 
              Table[
                Position[eq4, Derivative[i][y][x]], {i, 1, 
                  ord}]; 
t8 = 
              Position[Reverse[t7], {}]; 
If[
              Length[t8] > 0, 
              Table[p = Insert[p, 0, t8[[i]]], {i, 1, 
                  Length[t8]}]]; 
{ord, Drop[p, 1], 
              r} 
];

\(flst[lf_] := 
            Block[{fy, m}, 
m = Length[lf]; 
	
              fy = RotateLeft[lf, Floor[m/2]]; 

              RotateRight[Fourier[fy]/\@m, Floor[m/2]]];\)


coef[nn_] := 
          Module[{Snc, u, l, dd1, furelist, lk, fl1, furelist1, a1, a2, a3, 
              mn}, 
Snc[x_] = 
              If[Abs[x] ≤ 
(10^(-10)), 1, 
                Sin[x]/x]; 
u[k_, d_] = 
              If[k == 0, 0, 
                (E^(I π  k  d/(2.  nn + 1)))/
                  k]; 
l[d_] = 
              Table[u[k, d], {k, -nn, nn}]; 
dd1[
                x_] = π nn   Cos[
                  (1/2.)  (1 + nn) π x] Snc[((nn π x)/2.) ]/
                  Snc[((π x)/2.) ]; 
furelist[
                d_] := ((2.  nn + 1)/(2 I)) flst[l[d]] // 
                Chop; 
lk[d_] = 
              Table[((2.  k + d)/(2  nn + 1. )), {k, -nn, 
                  nn}]; 
fl1[d_] = 
              Map[dd1[#] &, lk[d]]; 
furelist1[d_] := 
              Map[dd1[#] &, lk[d]]; 
a1 = 
              Flatten[{-1, lk[-0.8], lk[-0.6], lk[-0.4], 
                  lk[-0.2], lk[0], lk[0.2], lk[0.4], lk[0.6], lk[0.8], 
                  1}]; 
a2 = 
              Flatten[{0, furelist[-0.8], furelist[-0.6], 
                  furelist[-0.4], furelist[-0.2], furelist[0], 
                  furelist[0.2], furelist[0.4], furelist[0.6], furelist[0.8], 
                  0}]; 
a3 = 
              Flatten[{dd1[-1], furelist1[-0.8], furelist1[-0.6], 
                  furelist1[-0.4], furelist1[-0.2], furelist1[0], 
                  furelist1[0.2], furelist1[0.4], furelist1[0.6], 
                  furelist1[0.8], dd1[1]}]; 
mn = 
              Re[Transpose[{a1, Transpose[{a2, a3}]}]] // 
                N; 
Interpolation[mn, 
              InterpolationOrder -> 12]
];

<< Utilities`FilterOptions`;

NLinearDSolve[eqns_, y_, {var_, a_, b_}, opt___]:=Module[{optNLinearDSolve,optlist,optErrorNLinearDSolve,irp,pi,pg,mxgp,sn,absrel, eq, ezr,sss,ord,p,f,k,kk,kkk,flist,fflist,ffflist,sh,nnc,rrnnc,tcunf,iop,ll,clist,k2,nn,l,rr,aux1,pop,aux2,mat,ss,solution,ssnew,eqss,res,sol,sollist,err,step,ldo,rrdo,aux1do,popdo,bnew,aux2do,matdo,ssdo,solutiondo,ssnewdo,eqssdo,resdo,soldo,sollistdo,errdo,norm,func,newfunc,wght,optfint,ezrnew,posit,posnew,ls,rs,yu1,yu1do,yu2,yu2do,yu3,yu3do,yu4,yu4do,l1,l1do,l2,l2do},

optfint = FilterOptions[FunctionInterpolation, opt];
optNLinearDSolve = {SymbolicSolution,  PrecisionInformation,  MaxGridPoints, SolutionNorm, 
                  Error, InterpolationOrder, InterpolationPrecision, 
                  AccuracyGoal, PrecisionGoal, InterpolationPoints, 
                  MaxRecursion, WorkingPrecision, Weight};
optlist = 
                Table[{opt}[[i]] // First, {i, 1, Length[{opt}]}];
optErrorNLinearDSolve = 
                Complement[optlist, optNLinearDSolve];
If[Length[optErrorNLinearDSolve] > 0, 
                Table[Message[NLinearDSolve::unknnopt, 
                    optErrorNLinearDSolve[[i]]], {i, 1, 
                    Length[optErrorNLinearDSolve]}]; Break[]];
ClearAttributes[NIntegrate, HoldAll];
{irp, pi, pg, mxgp, sn, absrel, 
                  wght} = {SymbolicSolution, PrecisionInformation, 
                      PrecisionGoal, MaxGridPoints, SolutionNorm, Error, 
                      Weight} /. {opt} /. Options[NLinearDSolve];
If[wght === Automatic, wght = 1];
If[(absrel === 
                      Relative)  ||  (absrel === Absolute), , 
                Message[ NLinearDSolve::errmes]; absrel = Relative];
       If[(irp === True) || (irp === 
                      False), , Message[NLinearDSolve::symsol]; 
                irp = False];	
If[(pi === True) || (pi === 
                      False), , Message[NLinearDSolve::precinf]; 
                pi = False];	
If[(sn === L2) || (sn === 
                      Uniform), ,  sn = L2; 
                Message[NLinearDSolve::nrmes]];
If[ pg  ===  Automatic,  
                pg  =  $MachinePrecision - 11, 
                If[Not[IntegerQ[pg]] || Not[Positive[pg]], 
                  Message[NLinearDSolve::prg]; pg = 5]];
If[ mxgp === Automatic, mxgp = 200, 
                If[Not[IntegerQ[mxgp]] || Not[Positive[mxgp]],  
                  mxgp  =  200; Message[NLinearDSolve::mxgpmes]]];
If[absrel === Relative, absrel = 1, 
                absrel = 2];
If[sn === L2, sn = 2, sn = 0];
pg = (10^(-pg));

 eq = eqns // First;   
 ezr = Drop[eqns, 1] // N;
sss = RecognizeDif[eq, y, var];
ord = sss[[1]];
If[Length[ezr] == ord, , 
                Message[NLinearDSolve::ndnco, Length[ezr], ord]; Break[]];
p[x_] = 
                sss[[2]] /. var -> x;
f[x_] = 
                sss[[3]] /. var -> x;
k[x_, 
                  t_] = (∑_(m = 1)^ord) p[
                            x][[m]]/((m - 
                              1)!) (((x - t))^(m - 1));
kk[x_, t_] = 
                k[(b - a) x + a, (b - a) t + a];
kkk = 
                Compile[{x, {t, _Real, 1}}, 
                  Evaluate[(b - a) 
                      kk[Sin[((π x)/2)]^2, 
                        Sin[((π t)/2)]^2]]];
flist[x_] = 
                Flatten[{f[x], 
                    Table[(∑_(s = m)^ord) p[
                              x][[s]] ((((x - a))^(s - m))/((s - m)!)), {m, 1, ord}]}];
fflist[x_] = 
                flist[(b - a) x + a];
ffflist[x_] = 
                fflist[Sin[((π x)/2)]^2];
SetAttributes[ffflist, Listable];
SetAttributes[fflist, Listable];
SetAttributes[flist, Listable];

sh = 110;

nnc = Round[(b - a) 31];
rrnnc = Range[1, nnc];
tcunf = 
                Table[ 2.  s/(2 nnc + 1), {s, 1, nnc}];
\(iop = \ 
                Sin[π\ \ tcunf]\ \ Apply[Plus, 
                    Map[\(\(1.  - \((\(-1\))\)\^#\)\/#\) 
                          Sin[π\ #\ \ tcunf] &, rrnnc]]\);
ll = 
                Table[(b - a) 
                      Sin[π/2  ((2.  s)/(2 nnc + 1))]^2 + a, {s, 
                    1, nnc}];
wlist = 
                Map[(wght /. var -> #) &, ll]; 
k2[x_, t_, 
                  s_] =  (((b - 
                          a))/((s - 
                            1)!)) (((x - (b - a) Sin[π/2 t]^2 - a))^(s - 1));

nn = 2;
l = 
                Table[((2.  s)/(2 nn + 1)), {s, 1, nn}];
rr = Range[1, nn];
aux1=Map[kkk[#, l]  Sin[π  l] &, 
                  l];
pop = 
                  Compile[{n, {x, _Real, 1}}, 
                    Evaluate[(1. /n) (1 - Cos[π n  l]) 
                        Sin[π n x]]];
aux2 = (2. /((2  nn + 1. ))) 
                    Apply[Plus, Transpose[Map[pop[#, l] &, rr]], {1}];
mat=IdentityMatrix[Length[l]] + aux1   aux2;
ss = 
                Transpose[LinearSolve[mat, ffflist[l]]];
ezr = 
                ezr /. Integrate -> NIntegrate;
pos = Position[ezr, NIntegrate];

If[
                Length[pos] > 0, 
ezrnew = 
                  Table[ezr[[i, 1]] - ezr[[i, 2]] + 0. , {i, 1, 
                      Length[ezr]}]; 
posit = 
                  Table[Position[ezrnew[[i]], y], {i, 1, 
                      Length[ezrnew]}]; 
posnew = 
                  Table[
                    Table[{posit[[j, i, 1]]}, {i, 1, 
                        Length[posit[[j]]]}], {j, 1, 
                      Length[ezrnew]}]; 
ls = -Apply[
                      Plus, Table[
                        Extract[ezrnew[[j]], posnew[[j]]], {j, 1, 
                          Length[ezrnew]}], 1]; 
rs = 
                  Table[Delete[ezrnew[[j]], posnew[[j]]], {j, 1, 
                      Length[ezrnew]}]
];

\(If[
                Length[pos] > 
                  0, 

clist = 
                  Table[c[m], {m, 1, ord}]; 
yu1[x_, 
                    s_] = \(2\/\((2\ nn + 1)\)\) 
                    ss . \((k2[x, l, s] Sin[π\ \ l] 
                          Apply[Plus, 
                            Map[\(1. \/#\) \((1 - 
Cos[π\ #\ \ ArcCos[1 - 
                                        2\ \(x - a\)\/\(b - a\)]/π])\) 
                                  Sin[π\ #\ \ l] &, 
                              rr]])\); 
yu2[x_, s_, m_] = 
                  Which[m == s, 1, m < s, \((x - a)\)\^\(s - m\), m > s, 0]\/\(\((s - m)\)!\); 
yu4[x_, 
                    s_] = \(yu1[x, 
                      ord - s]\)[\([1]\)]; 

yu3[x_, s_] = 
                  Table[yu2[x, ord - s, m], {m, 1, ord}] - 
                    Drop[yu1[x, ord - s], 1]; 
l1 = 
                  ls /. {\(\(Derivative[k_]\)[y]\)[cc_] -> yu4[cc, k], 
                      y[cc_] -> yu4[cc, 0]}; 
l2 = 
                  ls /. {\(\(Derivative[k_]\)[y]\)[cc_] -> yu3[cc, k], 
                      y[cc_] -> yu3[cc, 0]}; 
eqss = 
                  Table[\((l1 + l2 . clist)\)[\([i]\)] == 
                      rs[\([i]\)], {i, 1, ord}]; 
res = 
                  Solve[eqss, 
                    Intersection[
                      Extract[eqss, 
                        Position[eqss, c[i_]]]]]; 
ssnew[
                    x_, s_] = 
                  yu4[x, s] + clist . yu3[x, s], 
clist = 
                  Prepend[Table[\(-c[i]\), {i, 1, ord}], 
                    1]; 
\ \ \ solution[x_, 
                    s_] := 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2\/\((2\ nn + 1)\)\ clist . 
                        ss . \((k2[x, l, s] Sin[π\ \ l] 
                            Apply[Plus, 
                              Map[\(1. \/#\) \((1 - 
Cos[π\ #\ \ ArcCos[1 - 2\ \(x - a\)\/\(b - a\)]/π])\) 
                                    Sin[π\ #\ \ l] &, 
                                rr]])\) + ∑\+\(m = 1\)\%s\((\(c[
                                m]\/\(\((s - m)\)!\)\) 
                          If[m == s, 
                            1, \((x - a)\)\^\(s - m\)])\); 
ssnew[x_, s_] = 
                  solution[x, ord - s]; 
eqss = 
                  ezr /. {\(\(Derivative[k_]\)[y]\)[cc_] -> 
                        ssnew[cc, k], 
                      y[cc_] -> ssnew[cc, 0]}; 
res = 
                  Solve[eqss, 
                    Intersection[Extract[eqss, Position[eqss, c[i_]]]]]]\);

sol = 
                Compile[{{x, _Real, 1}}, 
                  Evaluate[ssnew[x, 0] /. res // First]];
sollist = sol[ll];
err = (10^11);
step = 1;
Do[Which[ 
           step == 1, 
                            nndo =  Which[   (10^10) <= 
                            err/pg,                     5 + Round[1.5  nn], 
                                                        (10^8) <= err/pg < (10^10),         5 +  
                            Round[1.4  nn], 
                                                        (10^6) <= err/pg < (10^8),           5 +  
                            Round[1.3  nn], 
                                                        (10^4) <= err/pg < (10^6),           5 + 
                            Round[1.2 nn], 
                                                        (10^2) <= err/pg < (10^4),           5 + 
                            Round[1.1  nn], 
                                                         err/pg < 
                            (10^2),                        nn + 6, 
                                                         err/pg < 
                            10 ,                         nn + 6 ], 
           step == 
                        2, 
                           nndo =  
                        Which[ 
                          (10^10) <= 
                            err/pg,                    2 + 
                            Round[1.04  nn], 
                                                    (10^8) <= err/pg < (10^10),         2 + 
                            Round[1.03  nn], 
                                                    (10^6) <= err/pg < (10^8),           2 + 
                            Round[1.02  nn], 
                                                    (10^4) <= err/pg < (10^6),           2 + 
                            Round[1.01  nn], 
                                                    (10^2) <= err/pg < (10^4),           nn + 
                            2, 
                                                     err/pg < (10^2),                        nn + 
                            2, 
                                                     err/pg < 10 ,                         nn + 
                            2 ]
];

If[nndo < mxgp - 0.4 mxgp, 
                      step = 1, step = 2];
ldo = 
                      Table[((2.  s)/(2 nndo + 1)), {s, 1, nndo}];
rrdo = Range[1, nndo];
aux1do=Map[kkk[#, ldo]  Sin[π  ldo] &, ldo];

If[nndo < sh,
                     popdo = 
                              Compile[{n, {x, _Real, 1}}, 
                                Evaluate[(1. /n) (1 - 
                                        Cos[π n  ldo]) 
                                    Sin[π n x]]];
                     aux2do = (2. /((2  nndo + 1. ))) 
                                Apply[Plus, 
                                  Transpose[
                                    Map[popdo[#, ldo] &, rrdo]], {1}]
,

                    bnew[x_, 
                                z_] = coef[nndo][
                                  Mod[z - x + 1, 2] - 1] + coef[nndo][
                                  Mod[z + x + 1, 2] - 1];
                    aux2do=Transpose[((1  )/((2  nndo + 
                                        1. )))Map[bnew[0, #]-bnew[1/π  ArcCos[1-2  ((Sin[ π/2  ldo])^2)],#]&,ldo]];

];
matdo=IdentityMatrix[Length[ldo]] + 
                          aux1do   aux2do;
ssdo = 
                      Transpose[LinearSolve[matdo, ffflist[ldo]]];
\(If[
                      Length[pos] > 
                        0, 
yu1do[x_, 
                          s_] = \(2. \/\((2\ nndo + 1)\)\) 
                          ssdo . \((k2[x, ldo, s] Sin[π\ \ ldo] 
                                Apply[Plus, 
                                  Map[\(1. \/#\) \((1 - 
Cos[π\ #\ \ ArcCos[1 - 
                                        2. \ \(x - a\)\/\(b - a\)]/π])\) 
                                        Sin[π\ #\ \ ldo] &, 
                                    rrdo]])\); 
yu2do[x_, 
                          s_, m_] = 
                        
                        Which[m == s, 1. , m < s, \((x - a)\)\^\(s - m\), m > s, 0]\/\(\((s - m)\)!\); 
yu4do[x_, 
                          s_] = \(yu1do[x, 
                            ord - s]\)[\([1]\)]; 
yu3do[x_, 
                          s_] = Table[yu2do[x, ord - s, m], {m, 1, ord}] - 
                          Drop[yu1do[x, ord - s], 
                            1]; 
l1do = 
                        ls /. {\(\(Derivative[k_]\)[y]\)[cc_] -> 
                              yu4do[cc, k], 
                            y[cc_] -> 
                              yu4do[cc, 0]}; 
l2do = 
                        ls /. {\(\(Derivative[k_]\)[y]\)[cc_] -> 
                              yu3do[cc, k], 
                            y[cc_] -> 
                              yu3do[cc, 0]}; 
eqssdo = 
                        Table[\((l1do + l2do . clist)\)[\([i]\)] == 
                            rs[\([i]\)], {i, 1, 
                            ord}]; 
resdo = 
                        Solve[eqssdo, 
                          Intersection[
                            Extract[eqssdo, 
                              Position[eqssdo, 
                                c[i_]]]]]; 
ssnewdo[x_, 
                          s_] = yu4do[x, s] + 
                          clist . 
                            yu3do[x, 
                              s], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ solutiondo[x_, 
                          s_] = 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(2\/\((2\ nndo + 1)\)\) 
                            clist . 
                              ssdo . \((k2[x, ldo, s] Sin[π\ \ ldo] 
                                  Apply[Plus, 
                                    Map[\(1. \/#\) \((1 - 
Cos[π\ #\ \ ArcCos[1 - 
                                        2\ \(x - a\)\/\(b - a\)]/π])\) 
                                        Sin[π\ #\ \ ldo] &, 
                                      rrdo]])\) + ∑\+\(m = 1\)\%s\((\(c[
                                      m]\/\(\((s - m)\)!\)\) 
                                If[m == s, 
                                  1, \((x - a)\)\^\(s - m\)])\); 

ssnewdo[x_, s_] = 
                        solutiondo[x, ord - s]; 
eqssdo = 
                        ezr /. {\(\(Derivative[k_]\)[y]\)[cc_] -> 
                              ssnewdo[cc, k], 
                            y[cc_] -> 
                              ssnewdo[cc, 
                                0]}; 

resdo = Solve[eqssdo, Intersection[
                            Extract[eqssdo, Position[eqssdo, c[i_]]]]]]\);


 soldo = 
                      Compile[{{x, _Real, 1}}, 
                        Evaluate[(ssnewdo[x, 0] /. resdo // First)]];
sollistdo = soldo[ll];
 If[
                      sn == 
                        2,   
            errdo = (((((2 (b - a))/(2 nnc + 1))))^(1/
                                2)) (((Apply[Plus, wlist  Abs[sollistdo - sollist]^2 iop]))^(1/2)); 
            norm = (((((2 (b - a))/(2 nnc + 1))))^(1/
                                2)) ((((Apply[Plus, wlist  Abs[sollistdo]^2  iop]))^(1/2))); 
            If[absrel == 1, 
                        errdo = 
                          errdo/norm ]
  , 
            errdo = 
                        Max[wlist  Abs[
                              sollistdo - 
                                sollist]]; 
            norm = Max[wlist  Abs[
                              sollistdo]]; 
            If[absrel == 1, 
                        errdo = errdo/norm ]
    ];
If[
                      errdo ≤ 3 pg, Break[]];
If[
                      nndo > mxgp, 
If[
                        absrel == 
                          1, 
                   If[sn == 2, 
                          Message[NLinearDSolve::fcrelL2, nn, norm, errdo, 
                            wght], 
                          Message[NLinearDSolve::fcrelun, nn, norm, errdo, 
                            wght]], 
                   If[sn == 2, 
                          Message[NLinearDSolve::fcabsL2, nn, norm, errdo, 
                            wght], 
                          Message[NLinearDSolve::fcabsun, nn, norm, errdo, 
                            wght]]]; 
Break[]];
nn = nndo;
sollist = sollistdo;
ssnew[x_, s_] = ssnewdo[x, s];

res = resdo;

ss = ssdo;
l = ldo;
rr = rrdo;
err = errdo, {50000000}];

 func[z_] = ssnew[z, 0] /. res // First;
 \(If[
                irp, 
\ \ \(Format[func] := 
                    SequenceForm["SymbolicSolution[", {a, 
                        b}, ",<>]"];\)\ \ \ 
, 
\(newfunc = 
                    FunctionInterpolation[Evaluate[func[x]], {x, a, b}, 
                      Evaluate[optfint], InterpolationOrder -> 6, 
                      InterpolationPrecision -> 
                        Min[pg + 1, $MachinePrecision], 
                      InterpolationPoints -> 200, 
                      PrecisionGoal -> Min[pg + 1, $MachinePrecision], 
                      AccuracyGoal -> 
                        Min[pg + 
                            1, $MachinePrecision]];\)\ \ \ \ \ \ \ \ \ \ 
]\);
 \(If[
                pi === True, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ If[sn == 
                    2, \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print["L2 norm (with weight->", 
                    wght, ") of solution -> about ", \ 
                    PaddedForm[norm, 
                      3]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ If[
                    absrel == 
                      1, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print["Relative L2  error (with weight->", 
                      wght, ")-> about ", 
                      PaddedForm[\ errdo, 
                        3]]\ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print["Absolute L2  error (with weight->", wght, ") -> about ", \ 
                      PaddedForm[\ errdo, 
                        3]]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print["Uniform norm (with weight->", 
                    wght, ") of solution -> about ", \ 
                    PaddedForm[norm, 
                      3]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ If[
                    absrel == 
                      1, \ 
\ \ \ \ \ \ \ \ \ \ \ \ Print["Relative Uniform  error (with weight->", wght, ") -> about ", 
                      PaddedForm[\ errdo, 
                        3]]\ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ \ Print["Absolute Uniform  error (with weight->", 
                      wght, ") -> about ", \ 
                      PaddedForm[\ errdo, 
                        3]]\ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]
\ \ \ ]\);
 If[
                irp,   Return[
                  func],  Return[newfunc];
 ];
];
posit = 
        Position[Options[NLinearDSolve], 
          SolutionNorm -> Uniform];
Options[NLinearDSolve] = 
        Drop[Options[NLinearDSolve], 
          posit // 
            First];
End[];
EndPackage[];

Chasing from MathSource  (Authors: Len Pismen, Leonid Braverman, Boris Rubinstein )

BeginPackage["Chasing`"]
Off[General::spell];
\(\(\(Off[General::spell1];\)\(
\)
    \)\)
\(\(\(AccuracyGoal::usage = "AccuracyGoal is an option for functions       Chasing corresponding to the       option AccuracyGoal of the function NDSolve.";\)\(
\)
    \)\)
\(\(\(PrecisionGoal::usage = "PrecisionGoal is an option for function       Chasing corresponding to the option PrecisionGoal of       the function NDSolve.";\)\(
\)
    \)\)
\(\(\(WorkingPrecision::usage = "WorkingPrecision is an option for function       Chasing corresponding to the option WorkingPrecision of       the function NDSolve.";\)\(
\)
    \)\)
\(\(MaxSteps::usage = "MaxStep is an option for function       Chasing corresponding to the option MaxSteps of       the function NDSolve.";\)

 (*Compiled::usage = "Compilied is an option for function         Chasing corresponding to the          option Compilied of the function NDSolve"; 
SolAccuracyGoal::usage = "SolAccuracyGoal         is an option for function Shooting corresponding to the          option AccuracyGoal of the function FindRoot."; 
SolWorkingPrecision::usage = "SolWorkingPrecision         is an option for function Shooting corresponding to the         option WorkingPrecision of the function FindRoot."; 
SolCompiled::usage = "SolCompilied is an option for function         Shooting corresponding to the         option Compilied of the function FindRoot.";*) 
\)
\(\(\(Method::usage = 
        Method::usage <> " When Method is used as an option for the function Chasing          it specifys the method of generation of the           auxillary problem. The default value Method -> Nonlinear          corresponds to the nonlinear problem.          The option value Method -> Linear corresponds to the linear          problem.";\)\(
\)
    \)\)
\(\(\(Linear::usage = "Linear is a value for option Method.";\)\(
\)
    \)\)
\(\(\(Nonlinear::usage = "Nonlinear is a value for option Method.";\)\(
\)
    \)\)
\(\(\(Chasing::usage = "Chasing[eqns,funcs,{x,xmin,xmax}]        finds a numerical solution to the linear boundary        value problem eqns for the functions funcs with the        independent variable x in the range xmin to xmax.       
Chasing[eqns,funcs,{x,xmin,xmax},{lam,lami,lamf}]        finds numerically an eigenvector of the       homogeneous linear boundary value problem eqns        for the functions funcs with the        independent variable x in the range xmin to xmax.       and the corresponding eigenvalue of       the parameter lam within the inteval {lami,lamf}.";\)\(
\)
    \)\)
\(\(\(Begin["`Private`"]\)\(
\)
    \)\)
{f, 
      g};   Off[General::spell, 
      General::spell1];
\(\(\(Chasing::pdif = "Equations `1` are        partial differential equations and cannot be        transformed.";\)\(
\)
    \)\)
\(\(\(Chasing::vform = "The supplied equations are not valid.";\)\(
\)
    \)\)
\(\(\(Chasing::vars = "Not all variables specified in the       list `1` appear in the equations.";\)\(
\)
    \)\)
\(\(\(Chasing::nonlin = "The specified problem is nonlinear. A linear       boundary value problem is expected.";\)\(
\)
    \)\)
\(\(\(Chasing::ndinf = "No infinite values allowed for the range.";\)\(
\)
    \)\)
\(\(\(Chasing::liminf = "No infinite values allowed in the boundary conditions.";\)\(
\)
    \)\)
\(\(\(Chasing::nend = "The backward integration cannot be finished. The result may       be inaccurate or cover only part of integration interval.";\)\(
\)
    \)\)
\(\(\(Chasing::eqnan = "Coefficients of differential equations       are not numbers.";\)\(
\)
    \)\)
\(\(\(Chasing::bcnan = "Coefficients of boundary conditions       are not numbers.";\)\(
\)
    \)\)
\(\(\(Chasing::nbcs = "Number of boundary conditions doesn't correspond to       differential order of equations or there are       boundary condition equations containing more than one point.";\)\(
\)
    \)\)
\(\(\(Chasing::ndnl = "Limits or boundary conditions points are       not real numbers.";\)\(
\)
    \)\)
\(\(\(Chasing::bmet = "Value of option Method -> `1` is not valid.";\)\(
\)
    \)\)
\(\(\(Chasing::endbc = "The result of integration may be inaccurate due       to possible incompatibility of boundary conditions.";\)\(
\)
    \)\)
\(\(\($Partder[
          equations_] := \((Message[Chasing::pdif, 
            equations]; {$Failed, $Failed, $Failed})\);\)\(
\)
    \)\)
\(\(\($NotValidForm := \((Message[
            Chasing::vform]; {$Failed, $Failed, $Failed})\);\)\(
\)
    \)\)
\(\(\($NotAllVars[
          vars_List] := \((Message[Chasing::vars, vars]; $Failed)\);\)\(
\)
    \)\)
\(\(\($InfLimits0 := \((Message[Chasing::ndinf]; $Failed)\);\)\(
\)
    \)\)
\(\(\($InfLimits := \((Message[Chasing::liminf]; $Failed)\);\)\(
\)
    \)\)
\(\(\($NonLinearProblem := \((Message[
            Chasing::nonlin]; $Failed)\);\)\(
\)
    \)\)
\(\(\($NoEndIntegration := Message[Chasing::nend];\)\(
\)
    \)\)
\(\(\($NonNumericEquation := \((Message[
            Chasing::eqnan]; $Failed)\);\)\(
\)
    \)\)
\(\(\($NonNumericBC := \((Message[Chasing::bcnan]; $Failed)\);\)\(
\)
    \)\)
\(\(\($BadBCNumber := \((Message[Chasing::nbcs]; $Failed)\);\)\(
\)
    \)\)
\(\(\($ComplexPts := \((Message[Chasing::ndnl]; $Failed)\);\)\(
\)
    \)\)
\(\(\($BadMethod[
          meth_] := \((Message[Chasing::bmet, meth]; $Failed)\);\)\(
\)
    \)\)
Clear[MinMaxDer];
\(\(\(MinMaxDer[vars_List, eqns_List] := 
      Module[{varpos, derpos, varsall, newvars1, newvarsall}, 
        varpos = 
          Map[Union[Cases[eqns, #[_] -> 0, Infinity]] &, 
            vars]; 
derpos = \((\(Cases[
                  eqns, \(\(Derivative[n_Integer]\)[#]\)[_] -> n, 
                  Infinity] &\) /@ vars)\); 
varsall = 
          Transpose[{vars, 
              MapThread[
                Join[#, #2] &, {varpos, 
                  derpos}]}]; 
{newvars1, newvarsall} = 
          Transpose[
            DeleteCases[
              varsall, {_, {}}]]; 
Transpose[{newvars1, 
            Map[Max, newvarsall]}]]\)\(
\)\(
\)
    \)\)
Clear[TimeDerivSelect];
\(\(\(TimeDerivSelect[eq_Equal, timevar_] := 
      Module[{eqq, timeord, timeordder}, 
        eqq = eq /. Equal -> Subtract; 
timeord = 
          Cases[Expand[eqq] + Null, _. \ \(\(Derivative[n__]\)[_]\)[args__] /; 
              Apply[Plus, 
                  Cases[Transpose[{{n}, {args}}], {l_Integer, 
                        timevar | timevar[_Integer]} -> l]] > 
                0]; 
timeordder = 
          Apply[Plus, timeord]; 
timeordder == 
          timeordder - eqq]\)\(
\)
    \)\)
\(\(TimeDerivSelect[eqs : {__Equal}, timevar_] := 
        Map[TimeDerivSelect[#, timevar] &, 
          eqs];\)

 (*Transformation\ of\ a\ system\ of\ ordinary\ differential\ equations\ of\ any\ order\ to\ a\ system\ of\ ODE' 
        s\ of\ a\ first\ order\ by\ introducing\ auxiliary\ dependent\ variables . 
          Boundary\ conditions\ are\ also\ \(\(transformed\)\(.\)\)*) 
\)
Off[Solve::svars];
Clear[TransformSystem];
\(\(\(TransformSystem[equations : {__Equal}, vars_List, t_] := 
      Module[{bc, eq, maxeq, newvars, lll, mapder, neweqs, rule0, rules, 
            eqs4, eqs, newbc, eqss, fs, rulefs, eqss1, eqss2}, 
          If[\(! FreeQ[
                equations, \(\(Derivative[__]\)[Alternatives @@ vars]\)[
                  B__ /; \((MemberQ[{B}, t] && Length[{B}] > 1)\)], 
                Infinity]\), 
            Return[$Partder[equations]]]; 
bc = 
            Cases[equations, a_ /; FreeQ[a, t]]; 
eq = 
            Complement[equations, 
              bc]; 
If[\(! FreeQ[
                eq, \((\((Alternatives @@ vars)\)[
                        aa_] | \(\(Derivative[_]\)[Alternatives @@ vars]\)[
                        aa_])\) /; UnsameQ[aa, t]]\), 
            Return[$NotValidForm]]; 
 (*Construction\ of\ the\ equations*) maxeq = MinMaxDer[vars, eq]; 
newvars = 
            MapThread[Table[\(#1[i]\)[t], {i, 0, #2 - 1}] &, 
              Transpose[maxeq]]; 
lll = 
            DeleteCases[
              Flatten[\(Partition[#, 2, 1] &\) /@ newvars, 
                1], {}]; 
mapder = \(MapAt[D[#, t] &, #, 
                  1] &\) /@ 
              lll; 
neweqs = \(\((# /. 
                    List -> Equal)\) &\) /@ 
              mapder; 
rule0 = \(\((#[aa_] -> \(#[0]\)[
                      aa])\) &\) /@ vars; 
rulemaxder = 
            Map[\((\(\(Derivative[Last[#]]\)[First[#]]\)[t] -> 
                    D[\(\(First[#]\)[Last[#] - 1]\)[t], t])\) &, 
              maxeq]; 
rules = \(\(\(Derivative[
                        nn_Integer]\)[#]\)[aa_] -> \(#[nn]\)[aa] &\) /@ 
              vars; 
eqs4 = \(\(eq /. rule0\) /. 
                rulemaxder\) /. rules; 
eqs = 
            TimeDerivSelect[Join[neweqs, eqs4], 
              t]; 
newbc = \(bc /. rule0\) /. 
              rules; 
eqss = 
            Join[eqs, newbc]; 
fs = 
            Array[f, 
              Apply[Plus, Map[Last, maxeq]]]; 
rulefs = 
            Thread[Map[Head, Flatten[newvars]] -> 
                fs]; 
eqss1 = 
            Map[Expand, eqs /. rulefs, {2}]; 
eqss2 = 
            Sort[eqss1, 
              First[Cases[#, \(\(Derivative[1]\)[_[n_Integer]]\)[_] -> 
                        n, Infinity, 1]] < 
                  First[Cases[#2, \(\(Derivative[1]\)[_[
                              n_Integer]]\)[_] -> n, Infinity, 
                      1]] &]; 
{Join[eqss2, 
              newbc /. rulefs], fs, rulefs}] /; 
        FreeQ[Names["System`*"], 
          Apply[Alternatives, Map[ToString, vars]]]\)\(
\)
    \)\)
Clear[SelectBC];
\(\(SelectBC[bc_List, boundPoints_List] := 
        Map[Cases[bc, a_ /; \(! FreeQ[a, _. \ \(_[_Integer]\)[#]]\)] &, 
          boundPoints];\)

\)
Clear[MakeAuxSystem]
\(\(\(MakeAuxSystem[actbc_, fs_List, funcs_List, matr_?MatrixQ, var_, 
        m_Integer] := 
      Module[{initc, pos, funcs1, add}, 
        initc = Join[
            Map[D[actbc /. \(f[n_]\)[_] -> f[n], #] &, 
              fs], {actbc /. \(f[n_]\)[_] -> 
                  0}]; 
pos = 
          Part[Position[Abs[Sign[initc]], 1, {1}, 1], 1, 
            1]; 
funcs1 = 
          ReplacePart[funcs, 
            m\ Part[initc, pos] + \((1 - m)\)\ Part[funcs, pos], 
            pos]; 
add = 
          m\ Part[matr, pos] . Drop[funcs, \(-1\)]/
                Part[initc, pos] /. {Part[funcs, pos] -> 
                1}; 
DeleteCases[
          Thread[\((D[funcs, var] + matr . Drop[funcs1, \(-1\)] - 
                  add\ funcs1)\) == 0], True]]\)\(
\)
    \)\)
Clear[Constructor]
\(\(\(Constructor[actbc_, fs_List, gs_List, var_] := 
      Module[{pt}, 
        pt = First[
            Cases[{actbc}, \(f[_]\)[a_] -> a, Infinity, 
              1]]; 
DeleteCases[
          Thread[\((gs /. var -> pt)\) == 
              Join[Map[D[actbc /. \(f[n_]\)[_] -> f[n], #] &, 
                  fs], {actbc /. \(f[n_]\)[_] -> 0}]], 
          True]]\)\(
\)\(
\)
    \)\)
Clear[EndPoint];
\(\(\(EndPoint[syst : {{__Equal} .. }, 
        funcs : {__List}, {var_, pt_, xmax_}, opts___?OptionQ] := 
      MapThread[
        EndPoint[#, #2, {var, pt, xmax}, opts] &, {syst, 
          funcs}]\)\(
\)\(
\)\( (*This\ old\ version\ is\ disabled\ EndPoint[syst : {__Equal}, funcs_?ListQ, {var_, pt_, xmax_}, 
            opts___?OptionQ] := 
        First[NDSolve[syst, funcs, {var, pt, xmax}, opts, 
              WorkingPrecision -> $MachinePrecision, 
              MaxSteps -> 2500, AccuracyGoal -> Automatic, 
              PrecisionGoal -> Automatic] /. 
            var -> 
              xmax]*) \)\(
\)\( (*This\ part\ is\ added\ 12\ Oct\ 1998. \ It\ solves\ the\ problem\ of\ divergence\ of\ the\ solution\ inside\ the\ interval\ \((pt, xmax)\) . 
          It\ checks\ whether\ the\ integration\ at\ this\ interval\ is\ finished\ successfully, 
      and\ if\ not\ it\ switches\ to\ integration\ along\ other\ contour\ \((function\ MakeCurve)\) . 
          This\ integration\ consists\ of\ three\ steps\ \(\((functions\ Step1, 2, 3)\)\(.\)\)*) \)\(
\)
    \)\)
\(\(\(EndPoint[syst : {__Equal}, funcs_?ListQ, {var_, pt_, xmax_}, 
        opts___?OptionQ] := 
      Module[{s0, rs}, 
        s0 = Step0[syst, funcs, {var, pt, xmax}, opts]; 
If[
          s0 === $Failed, 
          rs = Abs[xmax - pt]*0.1; 
MakeCurve[syst, 
            funcs, {var, pt, xmax}, rs, opts], 
          s0]]\)\(
\)\(
\)
    \)\)
\(\(\(Clear[Step0, Step1, Step2, Step3, MakeSyst, MakeCurve]\)\(
\)
    \)\)
\(\(\(Step0[syst : {__Equal}, funcs_?VectorQ, {var_, xmin_, xmax_}, 
        opts___?OptionQ] := 
      Module[{sol, endpt}, 
        Off[NDSolve::ndsz]; 
sol = 
          NDSolve[syst, funcs, {var, xmin, xmax}, opts, 
            WorkingPrecision -> $MachinePrecision, 
            MaxSteps -> 2500, AccuracyGoal -> Automatic, 
            PrecisionGoal -> Automatic]; 
endpt = 
          Part[sol, 1, 1, 2, 0, 1, 1, 2]; 
On[
          NDSolve::ndsz]; 
If[
          Chop[endpt - xmax] == 0, 
          First[sol /. 
              var -> 
                xmax], $Failed]]\)\(
\)\(
\)\( (*In\ order\ to\ make\ integration\ at\ steps\ 1 - 
          3\ we\ need\ to\ change\ the\ direction\ of\ integration\ at\ step\ 1\ x -> 
        xmin + I\ x; \((x = 
              0\ to\ rectside)\)\ at\ step\ 2\ x -> \(I\ rectside + 
            x\ \((x = xmin\ to\ xmax)\)\ at\ step\ 3\ x -> 
          xmax + I\ x\ \((x = 
                  rectside\ to\ 0)\)\)*) \)\(
\)
    \)\)
\(\(\(MakeSyst[syst : {__Equal}, funcs_?VectorQ, {var_, xmin_, xmax_}, 
        rectside_, step_Integer, prevres_List] := 
      Module[{eqs, bc, heads, rvar, lims, derfact, initpt, rin, rout}, 
        eqs = Select[syst, \(! FreeQ[#, var]\) &]; 
bc = 
          Complement[syst, eqs]; 
heads = 
          Map[Head, funcs]; 
{rvar, lims, derfact, 
            initpt} = 
          Switch[step, 
            1, {var -> xmin + I\ var, {var, 0, rectside}, \(-I\), 0}, 
            2, {var -> I\ rectside + var, {var, xmin, xmax}, 1, xmin}, 
            3, {var -> xmax + I\ var, {var, rectside, 0}, \(-I\), 
              rectside}]; 
rin = 
          Map[{#[var] -> #, \(#'\)[var] -> derfact*#'} &, heads] // 
            Flatten; 
rout = 
          Map[{# -> #[var], #' -> \(#'\)[var]} &, heads] // 
            Flatten; 
eqs = \(\(eqs /. rin\) /. rvar\) /. 
            rout; 
bc = 
          Map[\(Head[First[#]]\)[initpt] == Last[#] &, 
            prevres]; 
Join[eqs, bc]]\)\(
\)
    \)\)
\(\(\(Step1[syst : {__Equal}, funcs_?VectorQ, {var_, xmin_, xmax_}, 
        rectside_, opts___?OptionQ] := 
      Module[{syst1, sol, endpt}, 
        Off[NDSolve::ndsz]; 
syst1 = 
          MakeSyst[syst, funcs, {var, xmin, xmax}, rectside, 1, 
            Select[syst, FreeQ[#, var] &]]; 
sol = 
          NDSolve[syst1, funcs, {var, 0, rectside}, opts, 
            WorkingPrecision -> $MachinePrecision, 
            MaxSteps -> 2500, AccuracyGoal -> Automatic, 
            PrecisionGoal -> Automatic]; 
endpt = 
          Part[sol, 1, 1, 2, 0, 1, 1, 2]; 
On[
          NDSolve::ndsz]; 
If[
          Chop[endpt - rectside] == 0, 
          First[sol /. var -> rectside] // Chop, $Failed]]\)\(
\)
    \)\)
\(\(\(Step2[syst : {__Equal}, funcs_?VectorQ, {var_, xmin_, xmax_}, 
        rectside_, outstep1_, opts___?OptionQ] := 
      Module[{syst2, sol, endpt}, 
        If[outstep1 === $Failed, Return[$Failed]]; 
Off[
          NDSolve::ndsz]; 
syst2 = 
          MakeSyst[syst, funcs, {var, xmin, xmax}, rectside, 2, 
            outstep1]; 
sol = 
          NDSolve[syst2, funcs, {var, xmin, xmax}, opts, 
            WorkingPrecision -> $MachinePrecision, 
            MaxSteps -> 2500, AccuracyGoal -> Automatic, 
            PrecisionGoal -> Automatic]; 
endpt = 
          Part[sol, 1, 1, 2, 0, 1, 1, 2]; 
On[
          NDSolve::ndsz]; 
If[
          Chop[endpt - xmax] == 0, 
          First[sol /. var -> xmax] // Chop, $Failed]]\)\(
\)
    \)\)
\(\(\(Step3[syst : {__Equal}, funcs_?VectorQ, {var_, xmin_, xmax_}, 
        rectside_, outstep2_, opts___?OptionQ] := 
      Module[{syst3, sol, endpt}, 
        If[outstep2 === $Failed, Return[$Failed]]; 
Off[
          NDSolve::ndsz]; 
syst3 = 
          MakeSyst[syst, funcs, {var, xmin, xmax}, rectside, 3, 
            outstep2]; 
sol = 
          NDSolve[syst3, funcs, {var, rectside, 0}, opts, 
            WorkingPrecision -> $MachinePrecision, 
            MaxSteps -> 2500, AccuracyGoal -> Automatic, 
            PrecisionGoal -> Automatic]; 
endpt = 
          Part[sol, 1, 1, 2, 0, 1, 1, 1]; 
On[
          NDSolve::ndsz]; 
If[Chop[endpt] == 0, 
          Thread[\((Map[First, First[sol]] /. 
                    var -> xmax)\) -> \((Map[Last, First[sol]] /. 
                    var -> 0)\)] // Chop, $Failed]]\)\(
\)
    \)\)
\(\(\(MakeCurve[syst : {__Equal}, funcs_?VectorQ, {var_, xmin_, xmax_}, 
        rectside_, opts___?OptionQ] := 
      Module[{outstep1, outstep2}, 
        outstep1 = 
          Step1[syst, funcs, {var, xmin, xmax}, rectside, 
            opts]; 
outstep2 = 
          Step2[syst, funcs, {var, xmin, xmax}, rectside, outstep1, 
            opts]; 
Step3[syst, funcs, {var, xmin, xmax}, 
          rectside, outstep2, opts]]\)\(
\)
    \)\)
\(\(\(MakeMainSystem[groupauxsyst_List, groupgs_List, boundpts_List, 
        gs_List, lsortbounds_List, var_, xmax_, 
        newfuncs_List, {wp_, maxstep_, accgoal_, precgoal_, acc_}] := 
      Module[{finalgsrules, gmatr, mainmatr, mainfreeterm}, 
        finalgsrules = 
          Flatten[MapThread[
              EndPoint[#, #2, {var, #3, xmax}, WorkingPrecision -> wp, 
                  MaxSteps -> maxstep, AccuracyGoal -> accgoal, 
                  PrecisionGoal -> precgoal] &, {groupauxsyst, groupgs, 
                boundpts}]]; 
gmatr = \(gs /. 
              var -> xmax\) /. 
            finalgsrules; 
mainmatr = 
          Join[Map[Drop[#, \(-1\)] &, gmatr], 
            Outer[D, \(lsortbounds /. 
                  Equal -> Subtract\) /. \(f[n_]\)[_] -> f[n], 
              newfuncs]]; 
mainfreeterm = 
          Join[Map[Last, 
              gmatr], \(lsortbounds /. 
                Equal -> Subtract\) /. \(f[n_]\)[_] -> 
                0]; 
Chop[{mainmatr, mainfreeterm}, 
          acc]]\)\(
\)
    \)\)
Clear[Resolvent]
\(\(\(Resolvent[matr_?MatrixQ] := 
      Module[{es, ces, eval, efun, cefun}, 
        es = Chop[Eigensystem[matr]]; 
ces = 
          Chop[Eigensystem[Transpose[matr]]]; 
{eval, 
            efun} = Simplify[
            Transpose[
              DeleteCases[
                Transpose[es], {0, _List}]]]; 
cefun = 
          Simplify[
            Last[Transpose[
                DeleteCases[
                  Transpose[ces], {0, _List}]]]]; 
Apply[
          Plus, MapThread[\(Outer[Times, #, #2]/# . #2\)/#3 &, {efun, cefun, 
              eval}]]]\)\(
\)\(
\)
    \)\)
Clear[conjugate]
\(\(conjugate[expr_] := 
        expr /. Complex[a_, b_] -> 
            Complex[a, \(-b\)];\)

\)
Clear[Chasing];
\(\(\(Options[Chasing] = {Tolerance -> 10, 
          WorkingPrecision -> $MachinePrecision, MaxSteps -> 1000, 
          AccuracyGoal -> Automatic, PrecisionGoal -> Automatic, 
          Method -> Nonlinear};\)\(
\)
    \)\)
\(\(\(Chasing[eqns_List, 
        func : \((_Symbol | _Symbol[var_])\), {var_, xmin_?NumericQ, 
          xmax_?NumericQ}, opts___?OptionQ] := 
      Chasing[eqns, {func}, {var, xmin, xmax}, {0}, opts]\)\(
\)
    \)\)
\(\(\(Chasing[eqns_List, 
        funcs_List, {var_, xmin_?NumericQ, xmax_?NumericQ}, opts___?OptionQ] := 
      Chasing[eqns, funcs, {var, xmin, xmax}, {0}, opts]\)\(
\)
    \)\)
\(\(\(Chasing[eqns_List, 
        func : \((_Symbol | _Symbol[var_])\), {var_, xmin_?NumericQ, 
          xmax_?NumericQ}, lamset : {_Symbol, _?NumericQ, _?NumericQ}, 
        opts___?OptionQ] := 
      Chasing[eqns, {func}, {var, xmin, xmax}, {lamset}, 
        opts]\)\(
\)\(
\)
    \)\)
\(\(\(Chasing[eqns_List, 
        funcs_List, {var_, xmin_?NumericQ, xmax_?NumericQ}, 
        lamset : {_Symbol, _?NumericQ, _?NumericQ}, opts___?OptionQ] := 
      Chasing[eqns, funcs, {var, xmin, xmax}, {lamset}, 
        opts]\)\(
\)\(
\)
    \)\)
\(\(\(Chasing[eqns_List, 
        funcs_List, {var_, xmin_?NumericQ, xmax_?NumericQ}, 
        lamset : \(({{_Symbol, _?NumericQ, _?NumericQ} .. } | {0})\), 
        opts___?OptionQ] := 
      Module[{func1, eqns1, newfuncs, rulefs, needrulefs, dim, eqs, bounds, 
          boundpts, sortbounds, eqs0, testnum, eqnumtest, lamtestrule, 
          bcnumtest, boundlens, leadcoefs, jacob, jacobbc, tol, method, meth, 
          wp, maxstep, accgoal, precgoal, eigenflag, lamrule, lamset1, lam, 
          acc, freeterm, fullmatr, gs1, gs, auxsyst, groupauxsyst, groupgs, 
          finalgsrules, gmatr, mainmatr, mainfreeterm, endsol, newfs, 
          initcond, fsf, fsb, ic, fdrop, finalsol}, 
        func1 = Map[If[FreeQ[#, var], #, Head[#]] &, 
            funcs]; 
{eqns1, newfuncs, rulefs} = 
          TransformSystem[eqns, func1, var]; 
If[
          eqns1 === $Failed, 
          Return[$Failed]]; 
needrulefs = 
          Map[Reverse, 
            Cases[rulefs, Literal[_[0] -> _]]]; 
dim = 
          Length[newfuncs] + 1; 
eqs0 = 
          Cases[eqns1, 
            a_ /; \(! FreeQ[a, var]\)]; 
bounds = 
          Complement[eqns1, eqs0]; 
leadcoefs = 
          Map[First, eqs0] /. \(f[_Integer]'\)[var] -> 
              1; 
eqs = 
          Thread[Map[First, eqs0]/leadcoefs == 
              Map[Last, eqs0]/leadcoefs]; 
If[
          Abs[xmin] === Infinity || Abs[xmax] === Infinity, 
          Return[$InfLimits0]]; 
 (*We\ need\ check\ for\ numerical\ value\ of\ the\ equations\ coefficients*) testnum = 
          xmin + 0.1*\((xmax - xmin)\)*I; 
lamtestrule = 
          If[lamset =!= {0}, 
            Map[
              Part[#, 1] -> \((Part[#, 3] - Part[#, 2])\)*0.1*I + 
                    Part[#, 2] &, 
              lamset], {}]; 
eqnumtest = 
          N[\(\(Map[Last, eqs] /. \(f[_Integer]\)[var] -> 1\) /. 
                var -> testnum\) /. lamtestrule]; 
If[
          Not[Apply[And, Map[NumberQ, eqnumtest]]], 
          Return[$NonNumericEquation]]; 
boundpts = 
          Union[Cases[bounds, \(f[_Integer]\)[a_] -> a, 
              Infinity]]; 
If[
          Not[Apply[And, Map[NumberQ, N[boundpts]]]], 
          Return[$NonNumericBC]]; 
If[\(! FreeQ[
              N[Join[boundpts, {xmin, xmax}]], _Complex]\), 
          Return[$ComplexPts]]; 
xmaxbc = 
          Max[boundpts]; 
If[Abs[xmaxbc] === Infinity, 
          Return[$InfLimits]]; 
sortbounds = 
          SelectBC[bounds, boundpts]; 
bcnumtest = 
          N[\(\(Flatten[sortbounds] /. 
                  Equal -> Subtract\) /. \(f[_Integer]\)[_] -> 
                  testnum\) /. lamtestrule]; 
If[
          Not[Apply[And, Map[NumberQ, bcnumtest]]], 
          Return[$NonNumericBC]]; 
boundlens = 
          Map[Length, sortbounds]; 
If[
          Apply[Plus, boundlens] ≠ Length[eqs], 
          Return[$BadBCNumber]]; 
jacob = 
          Outer[D, Map[Last, eqs] /. \(f[n_Integer]\)[_] -> f[n], 
              newfuncs] // Transpose; 
jacobbc = 
          Outer[D, \(sortbounds /. 
                Equal -> Subtract\) /. \(f[n_Integer]\)[_] -> f[n], 
            newfuncs]; 
If[\(! FreeQ[jacob, 
                Alternatives @@ newfuncs]\) || \(! FreeQ[jacobbc, 
                Alternatives @@ newfuncs]\), 
          Return[$NonLinearProblem]]; 
method = \(\(Method /. {opts}\) /. Options[Chasing]\) /. {Linear -> 0, 
              Nonlinear -> 1}; 
If[
          method =!= 0 && method =!= 1, 
          Return[$BadMethod[
              method]]]; 
tol = \(Tolerance /. {opts}\) /. 
            Options[Chasing]; 
acc = 
          10^\((\(-tol\))\); 
{wp, maxstep, accgoal, 
            precgoal} = \({WorkingPrecision, MaxSteps, AccuracyGoal, 
                PrecisionGoal} /. {opts}\) /. 
            Options[Chasing]; 
freeterm = 
          Map[Last, eqs] /. \(f[n_Integer]\)[_] -> 
              0; 
eigenflag = 
          If[\((MatchQ[freeterm, {\((0)\) .. }])\) && lamset =!= {0}, 
            Length[lamset], 0]; 
fullmatr = 
          Join[jacob, {freeterm}]; 
gs1 = 
          Array[\(g[#, #2]\)[var] &, {dim - Last[boundlens] - 1, 
              dim}]; 
gs = 
          If[MatchQ[freeterm, {\((0)\) .. }] && 
              MatchQ[Map[Last, bounds], {\((0)\) .. }], 
            Map[ReplacePart[#, 0, \(-1\)] &, gs1], 
            gs1]; 
auxsyst = 
          MapThread[
              Join[MakeAuxSystem[#, newfuncs, #2, fullmatr, var, method], 
                  Constructor[#, newfuncs, #2, var]] &, {Flatten[
                    Drop[sortbounds, \(-1\)]] /. Equal -> Subtract, 
                gs}] // Chop; 
groupauxsyst = 
          First[Fold[{Append[First[#], Take[Last[#], #2]], 
                  Drop[Last[#], #2]} &, {{}, auxsyst}, 
              Drop[boundlens, \(-1\)]]]; 
groupgs = 
          First[Fold[{Append[First[#], Take[Last[#], #2]], 
                  Drop[Last[#], #2]} &, {{}, Map[DeleteCases[#, 0] &, gs]}, 
              Drop[boundlens, \(-1\)]]]; 
lamrule = 
          If[eigenflag == 0, {}, 
              ttts = Map[ToExpression[ToString[ttt] <> ToString[#]] &, 
                  Range[eigenflag]]; 
lamset1 = 
                MapThread[
                  ReplacePart[#, #2, 1] &, {lamset, 
                    ttts}]; 
lams = 
                Map[First, lamset]; 
Clear[
                foo]; 
foo[gh_] := 
                Module[{det}, 
                  det = Det[
                      First[\((Apply[
                                Function, {mmm[groupauxsyst, groupgs, 
                                      Drop[boundpts, \(-1\)], gs, 
                                      Last[sortbounds], var, xmaxbc, 
                                      newfuncs, {wp, maxstep, accgoal, 
                                        precgoal, acc}]} /. 
                                  Thread[
                                    lams -> Array[Slot, eigenflag]]] /. 
                              mmm -> MakeMainSystem)\) @@ 
                          gh]]; 
If[eigenflag == 1, 
                    Re[det], {Re[det], 
                      Im[det]}]]; 
\(FindRoot[foo[#], ##2, 
                    AccuracyGoal -> 15] &\) @@ Prepend[lamset1, ttts]] /. 
            Thread[ttts -> lams]; 
{mainmatr, 
            mainfreeterm} = 
          MakeMainSystem[groupauxsyst /. lamrule, groupgs, 
            Drop[boundpts, \(-1\)], gs, Last[sortbounds] /. lamrule, var, 
            xmaxbc, newfuncs, {wp, maxstep, accgoal, precgoal, 
              acc}]; 
endsol = 
          Which[MatchQ[mainfreeterm, {\((0)\) .. }] && 
              Chop[Det[mainmatr], acc] == 0, 
            Message[Chasing::endbc]; 
Chop[
              First[NullSpace[mainmatr, 
                  ZeroTest -> \((Chop[#, acc] == 0 &)\)]], acc], 
            MatchQ[mainfreeterm, {\((0)\) .. }] && 
              Chop[Det[mainmatr], acc] ≠ 0, 
            Message[Chasing::endbc]; 
mainfreeterm,  (*Not[
                  MatchQ[mainfreeterm, {\((0)\) .. }]] && 
                Chop[Det[mainmatr], acc] == 0, 
              Message[Chasing::endbc]; 
\(-Resolvent[
                    mainmatr] . mainfreeterm\), *) True, 
            Off[LinearSolve::nosol]; 
\(-LinearSolve[
                mainmatr, mainfreeterm]\)]; 
If[
          OffQ[LinearSolve::nosol], 
          On[LinearSolve::nosol]]; 
If[
          Head[endsol] =!= List, Return[{}]]; 
newfs = 
          Through[newfuncs[var]]; 
initcond = 
          Thread[\((newfs /. var -> xmaxbc)\) == 
              endsol]; 
 (*The\ backward\ integration\ process\ depends\ on\ the\ relation\ between\ xmax\ and\ xmaxbc . 
                There\ are\ three\ cases . 
                Case\ 1. \ \((xmax == xmaxbc)\) ==  > 
            backward\ integration\ is\ made\ from\ xmaxbc\ to\ xmin . 
                Case\ 2. \ \((xmax > xmaxbc)\) ==  > 
            Two\ function\ are\ produced - 
              one\ is\ by\ backward\ integration\ from\ xmaxbc\ to\ xmin, 
          and\ another\ by\ forward\ integration\ from\ xmaxbc\ to\ xmax . 
                Then\ these\ functions\ are\ merged . 
                Case\ 3. \ \((xmax < xmaxbc)\) ==  > 
            There\ are\ two\ backward\ integration\ \(steps : 
                first\ is\ from\ xmaxbc\ to\ xmax\); 
          the\ results\ are\ used\ as\ initial\ conditions\ for\ the\ second\ step - from\ xmax\ to\ \(\(xmin\)\(.\)\)*) finalsol = 
          Which[Chop[xmax - xmaxbc] == 0,  (*Case\ 1*) First[
              NDSolve[Join[eqs /. lamrule, initcond], 
                newfs, {var, xmax, xmin}, WorkingPrecision -> wp, 
                MaxSteps -> maxstep, AccuracyGoal -> accgoal, 
                PrecisionGoal -> precgoal]], 
            xmax > xmaxbc,  (*Case\ 2*) fsb = 
              First[NDSolve[Join[eqs /. lamrule, initcond], 
                  newfs, {var, xmaxbc, xmin}, WorkingPrecision -> wp, 
                  MaxSteps -> maxstep, AccuracyGoal -> accgoal, 
                  PrecisionGoal -> precgoal]]; 
fsf = 
              First[NDSolve[Join[eqs /. lamrule, initcond], 
                  newfs, {var, xmaxbc, xmax}, WorkingPrecision -> wp, 
                  MaxSteps -> maxstep, AccuracyGoal -> accgoal, 
                  PrecisionGoal -> 
                    precgoal]]; 
MapThread[
              MIF[#1, #2] &, {fsb, fsf}], 
            xmax < xmaxbc,  (*Case\ 3*) fsdrop = 
              First[NDSolve[Join[eqs /. lamrule, initcond], 
                  newfs, {var, xmaxbc, xmax}, WorkingPrecision -> wp, 
                  MaxSteps -> maxstep, AccuracyGoal -> accgoal, 
                  PrecisionGoal -> precgoal]]; 
If[
              Chop[Part[fsdrop, 1, 2, 0, 1, 1, 1] - xmax] ≠ 
                0,  (*Produce\ a\ warning\ message*) $NoEndIntegration; 
              Return[$Failed]]; 
ic = \(fsdrop /. 
                  var -> xmax\) /. 
                Rule -> Equal; 
First[
              NDSolve[Join[eqs /. lamrule, ic], newfs, {var, xmax, xmin}, 
                WorkingPrecision -> wp, MaxSteps -> maxstep, 
                AccuracyGoal -> accgoal, 
                PrecisionGoal -> precgoal]]];  (*hcihw*)  (*NOTE - 
            Backward\ integration\ may\ also\ fail\ due\ to\ sensitivity\ of\ the\ problem\ to\ the\ initial\ conditions, 
          divergency\ of\ the\ integration\ process\ itself, 
          and\ other\ reasons . 
              In\ order\ to\ be\ sure\ that\ the\ integration\ is\ finished, 
          we\ need\ to\ check\ whether\ xmin\ is\ reached*) If[
          finalsol === $Failed, Return[$Failed]]; 
If[
          Chop[Part[finalsol, 1, 2, 0, 1, 1, 1] - xmin] ≠ 0 || 
            Chop[Part[finalsol, 1, 2, 0, 1, 1, 2] - xmax] ≠ 
              0,  (*Produce\ a\ warning\ message*) $NoEndIntegration]; 
{Map[MakeF[#, Intersection[func1, funcs]] &, 
            Join[Select[finalsol /. needrulefs, FreeQ[#, f] &], 
                lamrule] /. \(f_[0]\)[a_] -> f[a]]}]\)\(
\)
    \)\)
Clear[MIF]
\(\(\(MIF[leftr_Rule, rightr_Rule] := 
      First[leftr] -> \(MergeIF[Part[leftr, 2, 0], Part[rightr, 2, 0]]\)[
          Part[leftr, 2, 1]]\)\(
\)
    \)\)
Clear[MergeIF]
\(\(\(MergeIF[left_InterpolatingFunction, right_InterpolatingFunction] := 
      InterpolatingFunction[{{left[\([1, 1, 1]\)], 
            right[\([1, 1, 2]\)]}}, {1, 
          Max[left[\([2, 2]\)], right[\([2, 2]\)]], 
          left[\([2, 3]\)] && 
            right[\([2, 3]\)], {False, Complex, 
              Real}[\([Min[{\(Position[{False, Complex, Real}, 
                      left[\([2, 4]\)]]\)[\([1, 
                      1]\)], \(Position[{False, Complex, Real}, 
                      right[\([2, 4]\)]]\)[\([1, 1]\)]}]]\)], 
          left[\([2, 5]\)], 
          left[\([2, 6]\)]}, {Join[left[\([3, 1]\)], 
            Drop[right[\([3, 1]\)], 1]]}, 
        Join[left[\([4]\)], 
          Drop[right[\([4]\)], 1]]]\)\(
\)\(
\)
    \)\)
Clear[MakeF]
\(\(\(MakeF[rl_Rule, heads_List] := 
      Map[Head, rl] /; MemberQ[heads, Head[First[rl]]]\)\(
\)
    \)\)
\(\(\(MakeF[rl_Rule, _List] := rl\)\(
\)
    \)\)
Clear[OffQ]
\(\(\(OffQ := \((Head[#] === $Off)\) &\)\(
\)\(
\)
    \)\)
End[];
EndPackage[];

Introduction

The package numerically solves a linear ordinary differential equation on a finite interval with conditions  given as linear combinations of the values of the solution and its derivatives in one or several points including integral distributions. The general problem has following form

(∑_(k = 0)^n)(Subscript[a, k])(x)(y^((k)))(x)=f (x);-∞<a≤x≤b<∞;(Subscript[a, n]) (x) ≠ 0, a<x<b ;                
  (∑_(p = 1)^((Subscript[m, s])))(Subscript[c, ps])(y^(((Subscript[k, ps]))))((Subscript[x, ps]))+(∑_(q = 1)^r)(∫_((Subscript[α, q]))^((Subscript[β, q])))(Subscript[z, q]) (t) y (t) \[DifferentialD]t= (Subscript[b, s]), s = 1,2,...,n ;
0 ≤ (Subscript[k, ps]) ≤ n - 1;  
a≤(Subscript[x, ps]),(Subscript[α, q]),(Subscript[β, q])≤b                                                                                                                                  (1)

with some constants {(Subscript[x, ps])},{(Subscript[c, ps])},{(Subscript[α, q])}, {(Subscript[β, q])}, {(Subscript[b, s])}, integers {(Subscript[k, ps])}, {(Subscript[m, s])},  r and functions {(Subscript[z, q]) (t)}, (Subscript[z, q]) (t) ⊂(Subscript[L, 1])((Subscript[α, q]),(Subscript[β, q])),when unique solvability is provided. Solution can be found with prescribed precision in a 
weighted (Subscript[L, 2]) (w (x), a, b)space.
    The name NLinearDSolve is conventional and may be changed. 
	The algorithm is best when both the coefficients {(Subscript[a, k])(x)}, functions {(Subscript[z, q]) (x)}, weight w(x) and the free term f(x) are  rather smooth and don't have sharp changing derivatives (i.e. have no high oscillations, peaks etc.). The program is not so sensitive to singularities on the end x=a and/or x=b if solution expected to be smooth (for example by physical point of view). In other cases the algorithm may works slow and sometimes cannot provide the prescribed precision.  
  The package works on the base of the numerical solution of the corresponding linear second kind integral equation on a finite interval which is essentially different from the source of current  MATHEMATICA  NDSolve package and the package Chasing of the MathSource. 
	Taking into account examples below we think that the best version for solution of linear ordinary differential equations may be the hybrid of  the packages NDSolve, Chasing and NLinearDSolve with their best properties
when linear equations are considered.

Implementation Notes

NLinearDSolve[eqns, y, {x, xmin, xmax}] finds a numerical solution to the spreaded given problems for ordinary linear differential equations eqns including pointwise or integral  conditions  for the function  y with the independent variable x in the range xmin to xmax. 
◼ NLinearDSolve  reduces initial, boundary and other spreaded given problems to the second kind Volterra type integral  equations. During recursion the number of interpolation  points  increases,  linear algebraic systems  solve and precision control works. 
◼ NLinearDSolve gives results in terms of  InterpolatingFunction or SymbolicSolution  objects. 
◼ If  SymbolicSolution is used as an output then NLinearDSolve gives values of the derivatives of the  solution in the ends of prscribed interval  using Limit function.
◼ SymbolicSolution  has a form of  a polynomial  which is represented using ChebyshevT polynomials. 
◼ SymbolicSolution objects are recommended to use when high precision is required. 
◼ The differential equations must contain enough initial, boundary or other conditions to determine the solutions correctly. 
◼ The differential equations in NLinearDSolve can involve complex numbers. 
◼  As a rule not more then prescribed precision is guaranteed for the solution but not for conditions.
◼ In addition to the options of FunctionInterpolation the following options can be given:

({{"\!\(\* StyleBox[\"\\\"SymbolicSolution\\\"\",\n\"MR\"]\) ", False"",            result is InterpolatingFunction  or SymbolicSolution object}, {"\!\(\* StyleBox[\"\\\"PrecisionInformation\\\"\",\n\"MR\"]\) ", "\!\(\*StyleBox[\"\\\"True\\\"\", \"MR\"]\) ", "           information about solution norm , discrepancy and error "}, {"\!\(\* StyleBox[\"\\\"MaxGridPoints\\\"\",\n\"MR\"]\) ", "\!\(\* StyleBox[\"\\\"Automatic\\\"\",\n\"MR\"]\) ", "           total number of grid "}, {"\!\(\* StyleBox[\"\\\"SolutionNorm\\\"\",\n\"MR\"]\) ", "\!\(\* StyleBox[\"\\\"Automatic\\\"\",\n\"MR\"]\)",            (Subscript[L, 2])  or   Uniform  norm and error are used}, {"\!\(\*
StyleBox[\"\\\"Error\\\"\",\n\"MR\"]\)
Weight               ", \(Automatic

          Automatic\ \*""   "\),           relative  or  absolute error to use
          weight in  (Subscript[L, 2]) space}})

◼ NLinearDSolve stops when  PrecisionGoal specified is met (default=$MachinePrecision-11). 
◼ The default setting of Automatic for MaxGridPoints is 200. 
◼ The default setting of  Automatic for SolutionNorm is (Subscript[L, 2]). 
◼ The default setting of  Automatic for Error is relative
◼ The default setting of  Weightt is w(x)≡1.

Demonstration of the package NLinearDSolve

Examples below are taken mainly from the well known books
1. E.KAMKE, Differentialgleichungen.Losunsmethoden und losunen, Leipzig , 1959.
2. H.Bateman,A.Erdelyi. Higher transcendental functions, vol.2,McGrow Hill,  1968

All examples are solved by the computer
 PENTIUM  4, 1.5 GHz, 512 RAM

1. Some comparisons with NDSolve and Chasing

1.1

Not all linear equations with  pointwise conditions can be solved by the NDSolve or Chasing .
 For example evaluatiion of following problems generates a lot of messages and no result.

NDSolve[{((y')')'[x]-x y'[x]+2 y[x]==-10 + x,y[-1] == 3,y'[2]==4,(y')'[4]==3},y,{x, -1, 5}]

Chasing[{((y')')'[x]-x y'[x]+2 y[x]==-10 + x,y[-1] == 3,y'[2]==4,(y')'[4]==3},y,{x, -1, 5}]

The same  problem can be easily solved by  NLinearDSolve

NLinearDSolve[{((y')')'[x]-x y'[x]+2 y[x]==-10 + x,y[-1] == 3,y'[2]==4,(y')'[4]==3},y,{x, -1, 5}]

Here is the plot of the solution

Plot[%[x], {x, -1, 5}, PlotRange -> All]

Non the less there are some examples when NDSolve or Chasing programs work faster than NLinearDSolve. Mainly it is connected with initial value problms.
 Here is one of the  simplest examples from the MATHEMATICA book. The evaluation of this problem  by NDSolve and Chasing works very fast

s1=NDSolve[{y'[x]==y[x],y[1] == 2},y,{x, 0, 1},PrecisionGoal -> 10,AccuracyGoal -> 12]//Timing

Plot[Abs[(y[x] /. s1[[2]]) - 2/E (E^x)], {x, 0, 1}, 
      PlotRange -> All]

s2=Chasing[{y'[x]==y[x],y[1] == 2},y,{x, 0, 1},PrecisionGoal -> 10,AccuracyGoal -> 12]//Timing

Plot[Abs[(y[x] /. s2[[2]]) - 2/E (E^x)], {x, 0, 1}]

Here NLinearDSolve works slower but a bit more precise

s3=NLinearDSolve[{y'[x]==y[x],y[1] == 2},y,{x, 0, 1},PrecisionGoal -> 10,SymbolicSolution -> True]//Timing

Plot[s3[[2]][x] - 2/E (E^x), {x, 0, 1}]

However if we are seeking the solution of the same problem on the interval (0, 5.5) instead of (0,1) then Chasing doesn't work

Chasing[{y'[x]==y[x],y[1] == 2},y,{x, 0, 5.5},PrecisionGoal -> 11,AccuracyGoal -> 11]

But our program works rather fast and generates demanded precision

s4=NLinearDSolve[{y'[x]==y[x],y[1] == 2},y,{x, 0, 5.5},PrecisionGoal -> 10,SymbolicSolution -> True,
PrecisionInformation -> True,Error -> Absolute]//Timing

Plot[s4[[2]][x] - 2 (E^(x - 1)), {x, 0, 5.5}, 
      PlotRange -> All]

Here NDSolve works very fast but can not reach to the prescribed precision

s5=NDSolve[{y'[x]==y[x],y[1] == 2},y,{x, 0, 5.5},PrecisionGoal -> 10,AccuracyGoal -> 15]//Timing

Plot[(y[x] /. s5[[2]]) - 2/E (E^x), {x, 0, 5.5}, 
      PlotRange -> All]

1.2

Not NDSolve nor Chasing can not solve previous equation with a simple boundary condition

NDSolve[{y'[x]==y[x],y[1] == y[0] + 1},y,{x, 0, 1}]

Chasing[{y'[x]==y[x],y[1] == y[0] + 1},y,{x, 0, 1}]

NLinearDSolve  generates the solution immediately and with high precision

s=NLinearDSolve[{y'[x]==y[x],y[1] == y[0] + 1},y,{x, 0, 1},SymbolicSolution -> True]//Timing

Plot[Abs[s[[2]][x] - ((E^x)/(E - 1))], {x, 0, 1}, 
      PlotRange -> All]

1.3

We are seeking a high-precision (with absolute uniform error approx. (10^(-12)) ) solution to the follwing boundary problem on the interval  (-1,3)

y''' (x)+ | x | 
        y (x) = (-1 + Abs[x]) BesselJ[
            1, (E^x)] - (E^(2 x)) (BesselJ[
                1, (E^x)] + 
              2 BesselJ[0, (E^x)] Sinh[
                  x]), 
y (-1) = 
      BesselJ[1, 
        (1/E)], 
y' (-1) = ((BesselJ[0, (1/E)] - BesselJ[2, (1/E)])/(2 E)), y (3) =   BesselJ[1, (E^3)]

with exact solution y(x)=BesselJ[1, (E^x)].

Here is the plot of the solution

Plot[BesselJ[1, (E^x)], {x, -1, 3}, 
      PlotRange -> All]

Below the program NLinearDSolve works much faster than  Chasing

w1 = 
      NLinearDSolve[{y'''[x] + 
                Abs[x] y[x] == (-1 + Abs[x]) BesselJ[
                    1, (E^x)] - (E^(2 x)) (BesselJ[1, (E^x)] + 
                      2 BesselJ[0, (E^x)] Sinh[x]), 
            y[-1] == 
              BesselJ[1, 
                (1/E)], y'[-1] == ((BesselJ[0, (1/E)] - BesselJ[2, (1/E)])/(2 E)), 
            y[3] ==  BesselJ[1, (E^3)]}, 
          y, {x, -1, 3}, 
PrecisionGoal -> 
            12, 
Error -> 
            Absolute, 
SolutionNorm -> 
            Uniform, 
SymbolicSolution -> True] // 
        Timing

Here is the shape of the SymbolicSolution

Take[w1[[2]][x] // Expand, 8]

Here is the plot of the absolute error

Plot[
        Evaluate[Abs[w1[[2]][x] - BesselJ[1, (E^x)]]], {x, -1, 3}, 
        PlotRange -> All] // Timing

Here are the discrepancies of the conditions

{w1[[2]][-1] - BesselJ[1, (1/E)], 
      Limit[D[w1[[2]][x], x], 
          x -> -1] - ((BesselJ[0, (1/E)] - BesselJ[2, (1/E)])/(2 E)), w1[[2]][3] -  
        BesselJ[1, (E^3)]}

The options in the program Chasing are chosen for achieving to the  required precision

w2 = 
      Chasing[{y'''[x] + 
                Abs[x] y[x] == (-1 + Abs[x]) BesselJ[
                    1, (E^x)] - (E^(2 x)) (BesselJ[1, (E^x)] + 
                      2 BesselJ[0, (E^x)] Sinh[x]), 
            y[-1] == 
              BesselJ[1, 
                (1/E)], y'[-1] == ((BesselJ[0, (1/E)] - BesselJ[2, (1/E)])/(2 E)), 
            y[3] ==  BesselJ[1, (E^3)]}, 
          y, {x, -1, 3}, PrecisionGoal -> 16, 
          AccuracyGoal -> 16, WorkingPrecision -> 22, 
          MaxSteps -> 2000] // Timing

Plot[
      Evaluate[Abs[(w2[[2]] // Last // Last // Last)[x] - 
            BesselJ[1, (E^x)]]], {x, -1, 3}, 
      PlotRange -> All]

This problem is unsolvable by NDSolve

w3 = 
      NDSolve[{y'''[x] + 
                Abs[x] y[x] == (-1 + Abs[x]) BesselJ[
                    1, (E^x)] - (E^(2 x)) (BesselJ[1, (E^x)] + 
                      2 BesselJ[0, (E^x)] Sinh[x]), 
            y[-1] == 
              BesselJ[1, 
                (1/E)], y'[-1] == ((BesselJ[0, (1/E)] - BesselJ[2, (1/E)])/(2 E)), 
            y[3] ==  BesselJ[1, (E^3)]}, 
          y, {x, -1, 3}] // Timing

2. Examples with pointwise given conditions that are unsolvable by  NDSolve and Chasing  as some  conditions are containing more than one point

2.1

Consider the following problem in the interval (2,3)

(x^2) y'' (x) + ((y (x))/(ln x)) = 
      x (e^x) (2 + x ln x), 
y (2) - 
        3  y' (2.3) + 1.5  y'' (3) = 
      1, 
y (2.1) + y' (3) = -1

Here is the solution

sol[x_] = (E^x) Log[x] - 19.864405477858284` Log[x] - 
          9.778619278759349` Log[
              x] (∫_2^x) 
                    (1/Log)[t]^2 \[DifferentialD]t;

NLinearDSolve solves this problem  rather fast and with high precision

w = 
      NLinearDSolve[{(x^2)  y''[x] + y[x]/Log[x] ==  
              x  (E^x) (2 + x  Log[x]), 
            y[2] - 3  y'[2.3] + 1.5  y[3] ==  
              1, 
y[2.1] + y'[3] == -1}, 
          y, {x, 2, 3}, 
PrecisionInformation -> 
            True, 
PrecisionGoal -> 
            14, 
SolutionNorm -> 
            L2, 
SymbolicSolution -> 
            True
] // Timing

Here are the actual (Subscript[L, 2])error and (Subscript[L, 2])norm of the solution which coincides with the information message

Sqrt[
      NIntegrate[Evaluate[Abs[w[[2]][x] - sol[x]]^2], {x, 2, 3}]/
        NIntegrate[Abs[sol[x]]^2, {x, 2, 3}]]

Sqrt[NIntegrate[Abs[sol[x]]^2, {x, 2, 3}]]

Here are the  discrepancies of the conditions

{y[2] - 3  y'[2.3] + 1.5  y[3] - 1, 
        y[2.1] + Limit[y'[x], x -> 3] + 1} /. 
      y -> w[[2]]

Evaluating this simple examples by NDSolve or  Chasing  generates  no result.

w = 
      NDSolve[{(x^2)  y''[x] + y[x]/Log[x] ==  
              x  (E^x) (2 + x  Log[x]), 
            y[2] - 3  y'[2.3] + 1.5  y[3] ==  
              1, 
y[2.1] + y'[3] == -1}, 
          y, {x, 2, 3}
] // Timing

w = 
      Chasing[{(x^2)  y''[x] + y[x]/Log[x] ==  
              x  (E^x) (2 + x  Log[x]), 
            y[2] - 3  y'[2.3] + 1.5  y[3] ==  
              1, 
y[2.1] + y'[3] == -1}, 
          y, {x, 2, 3}
] // Timing

2.2

Consider  the following numerically complicated problem in the interval ((10^(-7)) ,1) with a highly oscillating and non-smooth coefficient of the equation.

y'' (x)+ | sin ((1/x) - 1) | 
      y' (x) + (e^(2  x)) 
            y (x) == 
-28884 (E^(4 x)) cos (1 - 170 x) + (E^(6 x)) cos (1 - 
                170 x) + 
          1360 (E^(4 x)) sin (1 - 170 x) +  | 
      sin (1 - 
            (1/x)) |  (4 (E^(4 x)) cos (1 - 
                170 x) + 
          170 (E^(4 x)) sin (1 - 
                170 x)), 
y (0.0000001) + 
        3  y' (1) ==  
      17286.417630354244`, 
y (0.0000001) == 
      0.5403168269234947`

The solution of this problem also is  very oscillating function  y(x)=(E^(4 x))Cos[170 x - 1]   so the program works rather slow

w = 
      NLinearDSolve[{y''[x] + 
                Abs[Sin[1/ x - 1]] y'[x] + (E^(2  x)) 
                  y[x] == 
-28884 (E^(4 x)) Cos[1 - 170 x] + (E^(6 x)) Cos[
                    1 - 170 x] + 
                1360 (E^(4 x)) Sin[1 - 170 x] + 
                Abs[Sin[
                      1 - (1/x)]] (4 (E^(4 x)) Cos[
                          1 - 170 x] + 
                      170 (E^(4 x)) Sin[1 - 170 x]), 
            y[0.0000001] + 3  y'[1] ==  17286.417630354244`, 
            y[0.0000001] == 0.5403168269234947`}, 
          y, {x, 0.0000001, 1}, 
Error -> 
            Relative, 
PrecisionInformation -> 
            True, 
PrecisionGoal -> 
            11, 
SolutionNorm -> 
            Uniform
] // Timing

The result  is InterpolatingFunction object. Here is the graph of the relative  error. The precision of output is unsatisfactory as the function (E^(4 x)) Cos[170 x - 1]is highly oscillating and the values  of the options InterpolationPoints and InterpolationOrder are small.

Plot[
        Abs[w[[2]][x] - (E^(4  x)) Cos[170 x - 1]]/50.9, {x, 
          0.0000001, 1}, PlotRange -> All] // Timing

We increase the precision of the output by adding  the  options  InterpolationPoints->800 (default=200) and InterpolationOreder->13 (default=6)

w = 
      NLinearDSolve[{y''[x] + 
                Abs[Sin[1/ x - 1]] y'[x] + (E^(2  x)) 
                  y[x] == 
-28884 (E^(4 x)) Cos[1 - 170 x] + (E^(6 x)) Cos[
                    1 - 170 x] + 
                1360 (E^(4 x)) Sin[1 - 170 x] + 
                Abs[Sin[
                      1 - (1/x)]] (4 (E^(4 x)) Cos[
                          1 - 170 x] + 
                      170 (E^(4 x)) Sin[1 - 170 x]), 
            y[0.0000001] + 3  y'[1] ==  17286.417630354244`, 
            y[0.0000001] == 0.5403168269234947`}, 
          y, {x, 0.0000001, 1}, 
Error -> 
            Relative, 
PrecisionInformation -> 
            True, 
PrecisionGoal -> 
            11, 
SolutionNorm -> 
            Uniform, 
InterpolationPoints -> 
            800, 
InterpolationOrder -> 
            13
] // Timing

The result is good

Plot[
        Abs[w[[2]][x] - (E^(4  x)) Cos[170 x - 1]]/50.9, {x, 
          0.0000001, 1}, PlotRange -> All] // Timing

But the output can be given also as SymbolicSolution object which allows to avoid above mentioned problems much more fast

w = 
      NLinearDSolve[{y''[x] + 
                Abs[Sin[1/ x - 1]] y'[x] + (E^(2  x)) 
                  y[x] == 
-28884 (E^(4 x)) Cos[1 - 170 x] + (E^(6 x)) Cos[
                    1 - 170 x] + 
                1360 (E^(4 x)) Sin[1 - 170 x] + 
                Abs[Sin[
                      1 - (1/x)]] (4 (E^(4 x)) Cos[
                          1 - 170 x] + 
                      170 (E^(4 x)) Sin[1 - 170 x]), 
            y[0.0000001] + 3  y'[1] ==  17286.417630354244`, 
            y[0.0000001] == 0.5403168269234947`}, 
          y, {x, 0.0000001, 1}, 
Error -> 
            Relative, 
PrecisionInformation -> 
            True, 
PrecisionGoal -> 
            11, 
SolutionNorm -> 
            Uniform, 
SymbolicSolution -> 
            True
] // Timing

Plot[
        Evaluate[Abs[w[[2]][x] - (E^(4  x)) Cos[170 x - 1]]]/
          50.9, {x, 0.0000001, 1}, PlotRange -> All] // 
      Timing

2.3 Pointwise given problem for  8-th order differential equation

Consider the following  problem with spreaded pointwise conditionds on the interval (0,11)

y'''''''' (x) + x y' (x) - y (x) =  
      (1/128) (E^(sin (x))) (291 + 
            128 x cos (x) + 4192 cos (2 x) + 
            24220 cos (4 x) - 1056 cos (6 x) + cos (8 x) - 
            824 sin (x) - 19656 sin (3 x) + 8120 sin (5 x) - 
            56 sin (7 x)), 
y[0] == 
      sol[0], 
y[11] ==  
      sol[11], 
y[(11/2)] == 
      sol[(11/2)], 
y'[0] ==  sol'[
        0], 
y''[0] + 
        2 y[11/2] == sol''[0] + 
        2 sol[11/2], 
y'''[11] + y'[
          11/4] == sol'''[11] + sol'[
          11/4], 
y''''[0] == sol''''[
        0], y'''''[0] == sol'''''[0]

with  exact solution sol(x)=(e^(sin (x))).

sol[x_] = (E^Sin)[x]; a = 0; b = 11;

Here is the numerical solution

w=Timing[NLinearDSolve[{Derivative[8][y][x]+x y'[x]-y[x]==(1/128) (E^Sin)[x] (291 + 
                        128 x Cos[x] + 4192 Cos[2 x] + 24220 Cos[4 x] - 
                        1056 Cos[6 x] + Cos[8 x] - 824 Sin[x] - 
                        19656 Sin[3 x] + 8120 Sin[5 x] - 
                        56 Sin[7 x]),y[a] == sol[a],y[b] == sol[b],y[((a + b)/2)] == sol[((a + b)/2)],y'[a]==sol'[a],y''[a]+2 y[(b/2)]==sol''[a]+2 sol[(b/2)],Derivative[3][y][b]+y'[(b/4)]==Derivative[3][sol][b]+sol'[(b/4)],Derivative[4][y][a]==Derivative[4][sol][a],Derivative[5][y][a]==Derivative[5][sol][a]},y,{x, a, b}]]

Here is the actual (Subscript[L, 2])error

Sqrt[
      NIntegrate[Abs[w[[2]][x] - sol[x]]^2, {x, 0, 11}]/
        NIntegrate[Abs[sol[x]]^2, {x, 0, 11}]]

Here are the discrepancies of the conditions which are sensitive to derivatives used.  Note that the (Subscript[L, 2  ])-norm of the Derivative[5][y](x) is equal about 37 and the uniform error of it is equal about 26.

{y[a] - sol[a],y[b] - sol[b],y[((a + b)/2)] - sol[((a + b)/2)],y'[a]-sol'[a],y''[a]+2 y[(b/2)]-sol''[a]-2 sol[(b/2)],Derivative[3][y][b]+y'[(b/4)]-Derivative[3][sol][b]-sol'[(b/4)],Derivative[4][y][a]-Derivative[4][sol][a],Derivative[5][y][a]-Derivative[5][sol][a]}/.y -> w[[2]]

3. Boundary problem with complex numbers, weights and mixed pointwise and  integral type conditions

3.1

Consider the following problem on the interval (0,2)

Derivative[4][z][x]-Sin[x] Derivative[3][z][x]-I x z''[x]+Sin[x] z[x]=(E^(5 I x)) (-500 I + 635 x + 
              25 I (x^2) + (75 + (1 + 
                          125 I) x) Sin[x]),
z[(1/2)]-z'[(3/4)]=-1.723367395717695` + 
          3.947894980817675` I,
z[1]-3z'[0]=-3 + (E^(5 I)),
z'''[0] = -75,
(∫_0^2) z[
              x] \[DifferentialD]x = -0.29117370627792655` + 
          0.313873723763896` I

with  exact solution z[x] = x  (E^(5  I x)).

Here is the numerical solution

w=NLinearDSolve[{Derivative[4][y][x]-Sin[x] Derivative[3][y][x]-I x y''[x]+Sin[x] y[x]==(E^(5 I x)) (-500 I + 635 x + 
                        25 I (x^2) + (75 + (1 + 
                                    125 I) x) Sin[
                            x]),y[1]-3 y'[0]==-3 + (E^(5 I)),y[(1/2)]-y'[(3/4)]==-1.723367395717695 + 
                    3.947894980817675 I,Derivative[3][y][0]==-75,Integrate[
                    y[x], {x, 0, 2}] == -0.29117130551880593` + 
                    0.31386776719500614` I},y,{x, 0, 2},PrecisionInformation -> True]//Timing

Here is the absolute error

Plot[
      Evaluate[Abs[w[[2]][x] - x (E^(5  I x))]], {x, 0, 
        2}]

Here is the precision of the integral type condition. The actual  precision is much higher than it is required by the default settings.

Abs[
      NIntegrate[x (E^(5  I x)), {x, 0, 2}] - 
        NIntegrate[w[[2]][x], {x, 0, 2}]]

3.2

Consider  the following problem on the interval (1,2) with mixtured conditions

(x^(3 ))Derivative[4][y][x]+2 (x^2) Derivative[3][y][x]- x y''[x]+y'[x]-(x^3) y[x]==0,(∫_1^3) y[x] Sin[x] \[DifferentialD]x + y[1] == 
        9.42433203172718` + 
          4.65938872738068` I,(∫_1^2) y''[x] x \[DifferentialD]x + 
          2 (∫_1^3) y'[x] (x^2) \[DifferentialD]x == 
        77.16466705374886` + 
          41.744923777903146` I,y'[2] - 
          2 (∫_1^2) y''[x] 
                Log[x + 2] \[DifferentialD]x + (∫_2^3) y[
                x] \[DifferentialD]x == 
        2.962288718805786` + 
          2.377497782948788` I,y'''[1] == 
        1.7320907974056439` + 
          0.42941143422544675` I

Here is the exact solution

sol[x_] = 
        BesselJ[0, x] - BesselY[0, x] + 2 BesselJ[0, I x] + 
          BesselY[0, I x];

For solution of the problem we use both Integrate and NIntegrate functions. In the last one you can use the options of the NIntegrate function to calculate integrals more precise

w=NLinearDSolve[{(x^(3 ))Derivative[4][y][x]+2 (x^2) Derivative[3][y][x]- x y''[x]+y'[x]-(x^3) y[x]==0,(∫_1^3) y[x] 
                        Sin[x] \[DifferentialD]x + y[1] == 
                  9.42433203172718` + 
                    4.65938872738068` I,NIntegrate[y''[x] x, {x, 1, 2}, 
                      Method -> GaussKronrod, GaussPoints -> 11, 
                      PrecisionGoal -> 7] + 
                    2 (∫_1^3) y'[x] 
                          (x^2) \[DifferentialD]x == 
                  77.16466705374886` + 
                    41.744923777903146` I,y'[2] - 
                    2  NIntegrate[y''[x] Log[x + 2], {x, 1, 2}, 
                        Method -> GaussKronrod, GaussPoints -> 11, 
                        PrecisionGoal -> 7] + (∫_2^3) y[
                          x] \[DifferentialD]x == 
                  2.962288718805786` + 
                    2.377497782948788` I,y'''[1] == 
                  1.7320907974056439` + 
                    0.42941143422544675` I},y,{x, 1, 3},
PrecisionInformation -> True,
Weight -> (((x - 
                        2.5))^2) (((3 - x))^(1/2)),
SymbolicSolution -> True,
Error -> Absolute,
SolutionNorm -> Uniform,
PrecisionGoal -> 8]//Timing

Here is the graph of the absolute weighted error of the solution

Plot[
      Evaluate[(((x - 1.8))^2) (((3 - x))^(1/
                2)) (Abs[w[[2]][x] - sol[x]])], {x, 1, 3}, 
      PlotRange -> All]

Here is the graph the absolute not weighted error of the solution. The maximum of it  is about 2.5 times bigger.

Plot[Evaluate[(Abs[w[[2]][x] - sol[x]])], {x, 1, 3}, 
      PlotRange -> All]

You can use TrigExpand and Simplify functions to represent solution as an algebraic polynomial. This is helpful while calculating the derivatives of the solution

g[x_] = 
        Chop[w[[2]][x] // TrigExpand, (10^(-12))] // 
          Simplify;

g[x]

Here are the  discrepancies of the conditions

{NIntegrate[w[[2]][x] Sin[x], {x, 1, 3}] + w[[2]][1] - 
        9.42433203172718` -  4.65938872738068` I, 
      NIntegrate[w[[2]]''[x] x , {x, 1, 2}, 
          Method -> GaussKronrod, GaussPoints -> 11, 
          PrecisionGoal -> 7] + 
        2  NIntegrate[w[[2]]'[x] (x^2), {x, 1, 3}] - 
        77.16466705374886` - 
        41.744923777903146` I, g'[2] - 
        2  NIntegrate[g''[x] Log[x + 2], {x, 1, 2}, 
            Method -> GaussKronrod, GaussPoints -> 11, 
            PrecisionGoal -> 7] + 
        NIntegrate[w[[2]][x], {x, 2, 3}] - 2.962288718805786` -  
        2.377497782948788` I, g'''[1] - 
        1.7320907974056439` - 
        0.42941143422544675` I}

4. Singular Coefficients and general type conditions

4.1. Bessel's  equation

Consider the following problem with singular coefficients for well known Bessel equation on the interval (0,1)

(z^2) y''[z]+z y'[z]+((z^2) - (l^2)) y[z]=0,
 y (0) = 0, 
y (1) = 0.31166404072415316

with exact solution y(x)=BesselJ[1.3,z].

Here is the solution of the problem in the weighted space (Subscript[L, 2])((Sqrt[x]),0,1)

w=NLinearDSolve[{(z^2) y''[z]+z y'[z]+((z^2) - (1.3^2)) y[z]==0,y[0] == 0, y[1] == 0.31166404072415316},y,{z, 0, 1},
Weight -> (Sqrt[z]),
PrecisionInformation -> True]//Timing

Here is the Plot of the ggraph of the absolute error with the weight

Plot[(Sqrt[z]) Abs[w[[2]][z] - BesselJ[1.3, z]], {z, 0, 
        1}]

Here NDSolve and Chasing unsuccessfully attempt to solve this problem

NDSolve[{(z^2) y''[z]+z y'[z]+((z^2) - (1.3^2)) y[z]==0,y[0] == 0, y[1] == 0.31166404072415316},y,{z, 0, 1}]

Chasing[{(z^2) y''[z]+z y'[z]+((z^2) - (1.3^2)) y[z]==0,y[0] == 0, y[1] == 0.31166404072415316},y,{z, 0, 1}]

4.2. Equation for JacobiP[n,a,b,x] polynomials

Consider the following  problem on the interval (-1,1)

(1 - (x^2))w''[x]+(b - a - (a + b + 2) x)w'[x]+n (n + a + b + 1) w[x]=0,-1 ≤ x ≤ 
        1
w[0] - w'[0] + 
        2 (∫_(-1/2)^(1/2))w[x] 
              Sin[x] \[DifferentialD]x ==  
      0.19785845544020186`, 
w[-1] == -0.4656172500000001`, 
n = 
      5, 
a =  .2, b = - .3

with exact solution  w(x)=JacobiP[n,a,b,x].

Here is the  solution of the problem with the well known classic weight which is singular in the end x=-1 and wanishing in the end x=1

n = 5;
a =  .2;
b = - .3;
jac=NLinearDSolve[{(1 - (x^2))w''[x]+(b - a - (a + b + 2) x)w'[x]+
n (n + a + b + 1) w[x]== 0, w[-1] == -0.4656172500000001`,
w[0] - w'[0] + 
                    2  NIntegrate[Evaluate[w[x] Sin[x]], {x, 0, 1/2}, 
                        Method -> GaussKronrod, 
                        GaussPoints -> 
                          11] ==  -1.5239163778012081`},w,{x, -1, 1},
SymbolicSolution -> True,
Weight -> (((1 - 
                        x))^a) (((1 + x))^(b ))]//Timing

Here is the shape of the solution

jac[[2]][x] // Simplify

Here is a pure  polynomial form of the solution

j0 = jac[[2]][x] // TrigExpand

Here are the error function and its plot with the weight

j0 - JacobiP[n, a, b, x] // Expand

Plot[
      Evaluate[(((1 - x))^a) (((1 + x))^(b )) 
          Abs[jac[[2]][x] - 
              JacobiP[5,  .2, - .3, x]]], {x, -1, 1}, 
      PlotRange -> All]

4.3. The third order equation with singular coefficients

The general solution of the following considered equation has the form  (see below)

c[1] (E^x) + x c[2] BesselJ[1, I x] + 
      x c[3] BesselY[1, I x]

We are interesting for the solution

sol[x_] = 
      c[1] (E^x) + x c[2] BesselJ[1, I x] + 
          x c[3] BesselY[1, I x] /. {c[1] -> 0.5247162216356969`, 
          c[2] -> -0.746574013239847` + 
              1.2605358525241253` I, 
          c[3] -> -0.746574013239847` I}

Here NLinearDSolve works

w = 
      NLinearDSolve[{x  y'''[x] - 
                x y''[x] - (x - 1) y'[x] + (x - 1) 
                  y[x] == 0, 
y[0] ==  1, 
            y[1] == 1, y[2] == 0}, 
          y, {x, 0, 2}, 
PrecisionInformation -> 
            True, 
PrecisionGoal -> 
            4, 
SymbolicSolution -> 
            True
] // Timing

Here are the actual (Subscript[L, 2])error and (Subscript[L, 2])norm of the solution

Sqrt[
      NIntegrate[
        Evaluate[Abs[w[[2]][x] - sol[x]]^2], {x, 0, 2}]]

Sqrt[NIntegrate[Evaluate[Abs[sol[x]]^2], {x, 0, 2}]]

4.4. The same  equation with integral conditions

Let find the solution to the previous equation which is orthogonal to the subspace Span{Sin[π x],Cos[π x]}  and has  integral value equals to one. The exact
solution is

sol[x_] = 
        0.1684648737314839` (E^x) - (0.4904954997495638` - 
                0.26083376132874375` I) x BesselJ[
              1, I x] - 
          0.4904954997495643` I x BesselY[
              1, I x];

w = 
      NLinearDSolve[{x  y'''[x] - 
                x y''[x] - (x - 1) y'[x] + (x - 1) 
                  y[x] == 
              0, 
(∫_0^2) y[
                    x] \[DifferentialD]x ==  
              1, (∫_0^2) y[x] 
                  Sin[π x] \[DifferentialD]x == 
              0. , (∫_0^2) y[x] 
                  Cos[π x] \[DifferentialD]x == 0. }, 
          y, {x, 0, 2}, 
PrecisionInformation -> 
            True, 
PrecisionGoal -> 
            4, 
SymbolicSolution -> 
            True
] // Timing

Here is the actual (Subscript[L, 2])error and (Subscript[L, 2])norm of the solution

Sqrt[
      NIntegrate[
        Evaluate[Abs[w[[2]][x] - sol[x]]^2], {x, 0, 2}]]

Sqrt[NIntegrate[Abs[w[[2]][x]]^2, {x, 0, 2}]]

Here is discrepancies of the conditions

{NIntegrate[w[[2]][x], {x, 0, 2}] - 1, 
        NIntegrate[Sin[π x] w[[2]][x], {x, 0, 2}], 
        NIntegrate[Cos[π x] w[[2]][x], {x, 0, 2}]} // 
      N