Numerical Solution of 
One-Dimensional Linear Integral Equations of 
the Second Kind

Integral Equation Solver. Version Delta5

(*\(\(\(:\)\(Sources : 
	\ 1. \ L . M . Delves\)\)\  &\) J.L.Mohamed.
	  Computational Methods for Integral  Equations,
     Cambridge University Press,1985. 
  2.  A . Nersessian . 
                Generalization of the Euler-Maclaurin Formula and some Applications. Theory of functions and applications .  Collection of Works Dedicated to the Memory of 
Mkhitar M .  Djrbashian, "Louys" Publishing House, 1995, 133 - 138. 
	 *)
	  

	  

BeginPackage["NumericalMath`NISolve`"];
Off[
        General::spell1];
Off[General::spell];
Off[
        General::unfl];
Off[
        CompiledFunction::cfsa];
Off[
        CompiledFunction::cfn];
Off[
        Set::"setraw"];
Off[Set::write];
Off[
        Set::"wrsym"];
Off[
        FunctionInterpolation::"nreal"];
Off[
        NumberForm::"sigz"];
Off[
        ListInterpolation::"inddp"];
Off[
        General::"ovfl"];
Off[Break::"nofwd"];

NISolve::"usage" = "NISolve[eqn,sol,{x,xmin,xmax}] finds \
a numerical solution to the eqn for unknown function sol[x] with the \
independent variable x in the range xmin to xmax.";

NISolve::wneg = "Check your weight. It must be a non-negative \
function.";

NISolve::ncvb  = "NISolve failed to converge to prescribed \
precision after using `1` points of the grid. Current L_`2` norm (with weight \
`5`) of the solution equals about`3`, relative L_`2`-error (with weight `5`) \
equals about`4`. For more precision use more MaxGridPoints.";

NISolve::ncvb2  = "NISolve failed to converge to prescribed \
precision after using `1` points of the grid. Current L_`2` norm (with weight \
`5`) of the solution equals about`3`,absolute L_`2`-error (with weight `5`) \
equals about `4`. For more precision use more MaxGridPoints.";

NISolve:: 
            errmes = "Unknown option for Error. Now NISolve[] works with \
Error->Relative.";

NISolve::symsol = "Unknown option for SymbolicSolution. Now \
NISolve works with SymbolicSolution->False.";

NISolve::precinf = "Unknown option for PrecisionInformation. \
Now NISolve works with PrecisionInformation->False.";

NISolve::nrmes = "Unknown option for SolutionNorm. Now \
NISolve works with SolutionNorm->2.";

NISolve::comp = "Unknown option for Compiled. Now NISolve \
works with Compiled->True.";

NISolve::prg = "Unknown option for PrecisionGoal. Now NISolve \
works with PrecisionGoal->5.";

NISolve::mxgpmes = "Unknown option for MaxGridPoints. Now \
NISolve works with MaxGridPoints->200.";

NISolve::limitincons = "Your interval x∈[`1`,`2`] is \
inconsistent with the limits of integration.";

NISolve::snnr = "More than p=15 for L_p-norm is not \
recommended.";

NISolve::ncvb3  = "NISolve failed to converge to prescribed \
precision after using `1` points of the grid. Current L_`2` norm (with weight \
`5`) of the solution equals about`3`, relative L_`2`-error (with weight `5`) \
equals about`4`. For more precision use more MaxGridPoints or/and increase \
GaussPoints as well as MaxGridPoints.";

NISolve::ncvb4  = "NISolve failed to converge to prescribed \
precision after using `1` points of the grid. Current L_`2` norm (with weight \
`5`) of the solution equals about`3`,absolute L_`2`-error (with weight `5`) \
equals about`4`. For more precision use more MaxGridPoints or/and increase \
GaussPoints as well as MaxGridPoints.";

NISolve:: 
            mgp  = "More than 35 Gauss points of each panel are not \
recommended. Now NISolve[] works for GaussPoints->35.";

NISolve:: 
            lgp = "Less than 2 Gauss points on each panel are not \
recommended.  Now NISolve[] works for GaussPoints->2.";

NISolve:: 
            mgr =  "Your GaussPoints->`1` is more than MaxGridPoints->`2`. \
Now NISolve[] works for GaussPoints->`3`.";

NISolve::mprg = "More than 35 Gauss points of each panel are \
not recommended. Besides, your GaussPoints-> `1` is more than \
MaxGridPoints->`2`.  Now NISolve[] works for GaussPoints->`3`.";

NISolve::unknnopt1  = "Unknown Option `1`.";




  Options[NISolveFredholm] =Intersection[
            Flatten[{ {SymbolicSolution -> False,  Compiled -> True, 
                  GaussPoints -> Automatic, MaxGridPoints -> Automatic,  
                  PrecisionGoal -> Automatic,    PrecisionInformation -> 
                    False,  
SolutionNorm -> Automatic, 
                  Weight -> Automatic, Error -> Relative, 
                  Error -> Absolute}, Options[FunctionInterpolation]}]];

Options[NISolveVolterra]=Intersection[
            Flatten[{SymbolicSolution -> False, 
                Weight -> Automatic, PrecisionGoal -> Automatic, 
                MaxGridPoints -> Automatic, 
                PrecisionInformation -> False, 
                SolutionNorm -> Automatic, Error -> Relative, 
                Error -> Absolute, Options[FunctionInterpolation]}]];

Begin["NumericalMath`NISolve`Private`"];

<< Utilities`FilterOptions`;

NISolve[eq_, y_, {x_, c_, d_}, opt___] := 
          Module[{optFredholm, optlist, ss, k0, f0, optErrorFredholm, 
              optErrorVolterra}, 
optFredholm = {SymbolicSolution, Compiled, GaussPoints, MaxGridPoints,  
                PrecisionGoal,    PrecisionInformation, SolutionNorm, 
                Weight, Error, InterpolationOrder, InterpolationPrecision, 
                AccuracyGoal, InterpolationPoints, 
                MaxRecursion}; 
 
            optVolterra = {SymbolicSolution, Weight, PrecisionGoal, 
                MaxGridPoints, PrecisionInformation, SolutionNorm, Error, 
                InterpolationOrder, InterpolationPrecision, AccuracyGoal, 
                InterpolationPoints, 
                MaxRecursion}; 
optlist = 
              Table[{opt}[[i]] // First, {i, 1, 
                  Length[{opt}]}]; 
ss = 
              Recognize[eq, y, 
                x]; 
k0 = (Drop[ss, 1] // 
                      First)[[1]] // 
                First; 
f0 = (Drop[ss, 1] // First) // 
                Last; 
 
            If[ss[[1]] == 1, 
optErrorFredholm = 
                Complement[optlist, optFredholm]; 
If[
                Length[optErrorFredholm] > 0, 
                Table[Message[NISolve::unknnopt1, 
                    optErrorFredholm[[i]]], {i, 1, 
                    Length[optErrorFredholm]}]; 
                Return[$Failed]]; 
Return[
                NISolveFredholm[ k0,  
                  f0, {x, c, d},  {ss[[2, 2, 1]],  ss[[2, 3, 1]], 
                    ss[[2, 4, 1]]},  opt ]]]; 
If[
              ss[[1]] == 2, 
optErrorVolterra = 
                Complement[optlist, optVolterra]; 
If[
                Length[optErrorVolterra] > 0, 
                Table[Message[NISolve::unknnopt1, 
                    optErrorVolterra[[i]]], {i, 1, 
                    Length[optErrorVolterra]}]; 
                Return[$Failed]]; 
Return[
                NISolveVolterra[Drop[ss, 1] // First, {x, c, d}, 
                  opt]]];
];

\(Recognize[eq_, sol_, var0_] := 
	
          Module[{tt, toQ, lst, eqn, new, l, ll, ll0, tf, tk, lk, pos, 
              lev}, 
tt = Length[eq]; 
Do[
              eq = ReleaseHold[
                  eq /. {fac1_. \ \ Integrate[
                              int1_\ \ \ y[var1], {var1_, a1_, b1_}] + 
                          fac2_. \ \ Integrate[
                              int2_\ \ \ y[var2], {var2_, a1_, 
                                b1_}] -> 
HoldForm[
                          Integrate[\((fac1\ \ int1 + 
                                  fac2\ \ \((int2 /. var2 -> var1)\))\) 
                              y[var1], {var1, a1, 
                              b1}]]}], {tt}]; 

toQ[fac1_. \ Integrate[int1_, {var1_, a1_, b1_}], z_, 
                x_] := 
\ \ \ \ {Cancel[
                  fac1\ \ int1\ /\((z /. x -> var1)\)], {var1, a1, b1}}; 

            lst = Solve[eq, \ sol[var0]]; 
\ \ \ \ \ eqn = \(lst // First\) // 
                First; 
	\ \ \ new = {eqn[\([1]\)] == \((eqn[\([2]\)] // 
                        Expand)\)} // First; 
		l = 
              Position[new // Last, Integrate]; 
		ll = 
              Table[l[\([k]\)] // First, {k, 1, 
                  Length[l]}]; 
	\ \ \ \ \ \ \ \ ll0 = 
              Table[{l[\([k]\)] // First}, {k, 1, Length[l]}]; 
		tf = 
              Delete[new // Last, ll0]; 
\ 
            tk = If[Length[ll] == 1, {Part[\((new // Last)\), ll]}, 
                Part[\((new // Last)\), ll]]; 	
	
            lk = {Table[\(toQ[tk[\([i]\)], eqn[\([1]\)], 
                      var0]\)[\([1]\)], {i, 1, 
                    Length[tk]}], 
			Table[\(toQ[tk[\([i]\)], 
                      eqn[\([1]\)], var0]\)[\([2, 1]\)], {i, 1, 
                    Length[tk]}], 
				Table[\(toQ[tk[\([i]\)], 
                      eqn[\([1]\)], var0]\)[\([2, 2]\)], {i, 1, 
                    Length[tk]}], 
				Table[\(toQ[tk[\([i]\)], 
                      eqn[\([1]\)], var0]\)[\([2, 3]\)], {i, 1, 
                    Length[tk]}]}; 

            pos = Append[lk, tf]; 
If[
              Length[pos[\([1]\)]] > 1, 
              lev = 2, 
If[\((NumberQ[\(pos[\([3]\)]\)[\([1]\)] // 
                          N] || \((\(pos[\([3]\)]\)[\([1]\)] // 
                            N)\) === ∞ || \((\(pos[\([3]\)]\)[\([1]\)] // N)\) === \(-∞\))\) && \((NumberQ[\(pos[\([4]\)]\)[\([1]\)] // 
                          N] || \((\ \(pos[\([4]\)]\)[\([1]\)] // 
                            N)\) === ∞ || \((\ \(pos[\([4]\)]\)[\([1]\)] // N)\) === \(-∞\))\), lev = 1, 
                lev = 2]]; 

{lev, 
              pos}
];\)

\(NISolveFredholm[\ k0_, \ 
            f0_, {x_Symbol, c_, d_}, \ {y_Symbol, \ a0_, \ b0_}, \ 
            opt___\ ] := 
	
          Module[{fi, ww, pg, ss, sol, fintopt}, 
fintopt = 
              FilterOptions[FunctionInterpolation, 
                opt]; 
{fi, ww, 
                pg} = \({SymbolicSolution, Weight, PrecisionGoal} /. {opt}\) /. 
                Options[NISolveFredholm]; 
exmp = 
              If[ww === Automatic, 1, 
                ww]; 
If[\((fi === True)\) || \((fi === 
                    False)\), , Message[NISolve::symsol]; 
              fi = False]; 
\ \ \ If[\ 
              pg\  === \ Automatic, \ pg\  = \ $MachinePrecision - 11, 
              If[Not[IntegerQ[pg]] || Not[Positive[pg]], 
                Message[
                  NISolve::prg]]]; 

If[\((a0 === \(-Infinity\))\) || \((b0 === 
                    Infinity)\) || \((c === \(-Infinity\))\) || \((d === 
                    Infinity)\), 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ 	\ \ \ \ \ \ a = ArcTan[a0] // N; 
			\ \ b = 
                ArcTan[b0] // 
                  N; 
			\ \ k = \((k0 /. {x -> Tan[x], y -> Tan[y]})\)/
                  Cos[y]^2; 
			\ \ f = \((f0 /. 
                    x -> Tan[x])\); 
\ \ \ \ \ \ \ \ \ \ \ \ \ If[
                ww === Automatic, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ sol = 
                  NISolveFredholmLimit[\ k, \ 
                    f, \ {x, ArcTan[c], ArcTan[d]}, \ {y, \ a, \ b}, \ opt, 
                    Weight -> 
                      1/Cos[x]^2]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ If[fi, \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ss[
                      z_] = sol /. 
                      x -> 
                        ArcTan[
                          z]; \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Format[ss] := 
                    SequenceForm["SymbolicSolution[", {c, 
                        d}, ",<>]"]; \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                    
                    ss]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ sol = 
                    FunctionInterpolation[
                      Evaluate[sol], {x, ArcTan[c], ArcTan[d]}, 
                      Evaluate[
                        fintopt], 
\ \ \ \ \ \ \ \ \ \ \ InterpolationPrecision -> 
                        Min[pg + 
                            1, $MachinePrecision], 
\ \ \ \ \ \ \ \ \ \ \ \ InterpolationOrder -> 7, 
                      PrecisionGoal -> Min[pg + 1, $MachinePrecision], 
                      AccuracyGoal -> 
                        Min[$MachinePrecision, 
                          pg + 1], \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ InterpolationPoints -> 200, 
                      MaxRecursion -> 
                        50]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ss[z_] = sol[ArcTan[
                        z]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Format[ss] := SequenceForm["InterpolatingFunction[", {{c, 
                          d}}, ",<>]"]; \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                    ss]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]
\ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ sol = 
                  NISolveFredholmLimit[\ k, \ 
                    f, \ {x, ArcTan[c], ArcTan[d]}, \ {y, \ a, \ b}, \ 
                    Weight -> \((ww /. x -> Tan[x])\) 1/Cos[x]^2, 
                    opt]; 
		\ \ \ \ \ \ \ If[
                  fi, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ss[z_] = sol /. x -> 
                        ArcTan[
                          z]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Format[ss] := SequenceForm["SymbolicSolution[", {c, 
                        d}, ",<>]"]; \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                    ss]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ sol = 
                    FunctionInterpolation[
                      Evaluate[sol], {x, ArcTan[c], ArcTan[d]}, 
                      Evaluate[\ fintopt], 
                      InterpolationPrecision -> 
                        Min[pg + 
                            1, $MachinePrecision], 
\ \ \ \ \ \ \ \ \ \ \ \ InterpolationOrder -> 7, 
                      PrecisionGoal -> Min[pg + 1, $MachinePrecision], 
                      AccuracyGoal -> 
                        Min[$MachinePrecision, 
                          pg + 1], \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ InterpolationPoints -> 200, 
                      MaxRecursion -> 
                        50]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ss[z_] = 
                    sol[ArcTan[
                        z]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Format[ss] := SequenceForm["InterpolatingFunction[", {{c, 
                          d}}, ",<>]"]; \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                    ss]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\ 
\ \ \ \ \ \ ]
, 		
\ \ \ \ \ \ 
\ \ sol = 
                NISolveFredholmLimit[\ k0, \ 
                  f0, \ {x, c, d}, \ {y, \ a0, \ b0}, \ 
                  opt\ ]; 
\ \ \ \ \ \ \ \ If[
                fi, \ 
\ \ \ \ \ \ \ \ \ \ ss[z_] = 
                  sol /. x -> 
                      z; 
\ \ \ \ \ \ \ \ \ \ Format[ss] := 
                  SequenceForm["SymbolicSolution[", {c, 
                      d}, ",<>]"]; \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ Return[
                  ss]
\ , 
\ \ \ \ \ \ \ \ \ sol = FunctionInterpolation[Evaluate[sol], {x, c, d}, 
                    Evaluate[fintopt], \ \ \ \ InterpolationPrecision -> 
                      Min[pg + 
                          1, $MachinePrecision], 
\ \ \ \ \ \ \ \ \ \ \ \ InterpolationOrder -> 7, 
                    PrecisionGoal -> Min[pg + 1, $MachinePrecision], 
                    AccuracyGoal -> 
                      Min[$MachinePrecision, 
                        pg + 1], \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ InterpolationPoints -> 200, 
                    MaxRecursion -> 
                      50]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                  sol]
\ \ \ \ \ \ \ ]\ \ \ \ \ \ \ \ \ \ 
]

		];\)


NISolveFredholmLimit[ k0_,  
          f0_,  {x_Symbol, ccc_, ddd_},  {y_Symbol,  a0_,  b0_},  
          opt___ ] :=

   
 
   Block[\({\ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ Local\ variables\ *) \ 
 (*\ 
                For\ options\ *) 
	\ \ \ \ \ \ ag, pi, pg, sn, sn0, wp, 
              mxgp, n, 
              absrel, 

	 (*For\ the\ kernel, \ 
                free\ term\ and\ range*) 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ k, f, kc, fc, comp, a, b, aa, bb, ap, 
              aq, 
	
	 (*\ 
                Connected\ with\ Gauss\ quadratures\ *) 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ t, t0, jt, jtc, jtl, jtr, w, w0, jw, jwc, 
              jwl, jwr, 
	
	 (*Connected\ with\ LinearSolve[]*) 
	\ \ \ \ \ \ \ \ \ \ \ \ \ gh, sh, la, mmc, mmc1, 
              s, 
	
	 (*For\ morms\ of\ solution\ and\ error\ estimations\ *) 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ppr, eer, dd, phi, 
              wt, 
	
	 (*For\ following\ automatic\ chosing\ of\ final\ panels*) 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ u, v, 
              uu, uv, 
	
	 (*\ 
                Remaining\ *) 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ e, 
              off1, off2, off3, OnQ, ier, co, hi, hi1, q
\ \ \ \ }\),
		
		  

		OnQ[mess_] := 
                MemberQ[{MessageName,  String},  Head[mess]];
		off1 = OnQ[LinearSolve::luc];
		     off2 = OnQ[CompiledFunction::cfn];
		         off3 = OnQ[General::spell];
                           off4 = OnQ[NumberForm::sigz];
                                  off5 = 
                OnQ[NIntegrate::inum];
                                        off6 = 
                OnQ[Set::write];
                                              off7 = 
                OnQ[NIntegrate::nlim];

		If[off1, Off[LinearSolve::luc]];
		    If[off2, Off[CompiledFunction::cfn]];
                          If[off3, 
                Off[General::spell]];
                         If[off4, 
                Off[NumberForm::sigz]];
                              If[off5, 
                Off[NIntegrate::inum]];
                                   If[off6, 
                Off[Set::write]];
                                        If[off7, 
                Off[NIntegrate::nlim]];

ag = 1;
	         kk := Function[{x, y}, k0];
      ff := Function[{x}, f0];
    
	\({fi, pi, \ pg, sn, \ \ mxgp, \ comp, \ n, wght, 
                  absrel}\  = 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \({SymbolicSolution, PrecisionInformation, PrecisionGoal, 
                      SolutionNorm, \ \ MaxGridPoints, Compiled, GaussPoints, 
                      Weight, Error}\  /. \ {opt}\)\  /. 
                  Options[NISolveFredholm]\);
  
   
	

     

	   If[(absrel === Relative)  ||  (absrel === 
                      Absolute), , Message[ NISolve::errmes]; 
                absrel = Relative];
  	 If[(fi === True) || (fi === False), , 
                Message[NISolve::symsol]; fi = False];	
     If[(pi === True) || (pi === False), , 
                Message[NISolve::precinf]; pi = False];	
     If[wght === Automatic, wght = 1];
     If[sn ===  Automatic,  sn = 2];
     If[(sn === ∞) || (sn === 
                      Automatic), , 
                If[Not[IntegerQ[sn]] || Not[Positive[sn]], 
                  Message[NISolve::nrmes]; sn = 2]];
              	Which[sn < ∞, sn0 = sn;  
                If[sn > 15, Message[NISolve::snnr]], sn === Infinity,  
                sn0 = Infinity];
           If[(comp === True) || (comp === 
                      False), , Message[NISolve::comp]; comp = True];	


      mxgp  =  Max[ mxgp,  3];

     If[ (n  ===  Automatic)  ||  ! IntegerQ[
                      n], n  = 5 ];
     If[ pg  ===  Automatic,  pg  =  5, 
                If[Not[IntegerQ[pg]] || Not[Positive[pg]], 
                  Message[NISolve::prg]; pg = 5]];
     If[ (! IntegerQ[mxgp]) || (mxgp  ===  
                      Automatic),  mxgp  =  n  Round[200/n]];
      {ag,  pg}  =  10^(-{ag,  pg});
					
      If[
                n < 2
, 
n = 2;  
                Message[ NISolve::lgp]];
                  \(If[\ 
                n <= 35
\ \ \ \ \ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ If[mxgp < n, 
                  Message[NISolve::mgr, \ \ \ n, mxgp, mxgp]; \ \ \ \ \ n = 
                    mxgp; mxgp = mxgp + 1, True], 
	\ \ \ \ \ \ \ \ If[
                  mxgp > n, Message[NISolve::mgp]
	, 

                  Message[NISolve::mprg, n, mxgp, 
                    Min[mxgp, 35]]; 
	\ \ \ n = Min[mxgp, 35]; 
                  mxgp = mxgp + 1]]\);

 comp  =  If[ comp,  Compile,  Function ];
   
          a := a0;b := b0;
         kc    =  
                comp[ {x,  y}, Evaluate[ k0 ]];
         fc    =  comp[{x}, Evaluate[ f0]];
     
      
     {it, iw} = 
                Transpose[
                  Re[ GQW[ n, 0. ,  1. ,  $MachinePrecision ] ]];
      np = Round[mxgp/n];
      l =  -Log[10, pg] // N;
      ppr =  .85  pg;
     c0 = Ceiling[91/n];
     jtc = 
                Flatten[
                  Table[ccc + (ddd - ccc) (it + s)/c0, {s, 1, 
                      c0 - 2}]];
     jwc = 
                Flatten[Table[iw (ddd - ccc)/c0, {c0 - 2}]]; 
	  q = Ceiling[11/n] + 3;
 	ap = ccc + (ddd - ccc)/c0;
	 jtl = 
                Flatten[
                  Table[ccc + (ap - ccc) (it + s)/q, {s, 0, 
                      q - 1}]];
     jwl = Flatten[Table[iw (ap - ccc)/q, {q}]]; 
	aq = ccc + (ddd - ccc)/c0 (c0 - 1);
	jtr = 
                Flatten[
                  Table[aq + (ddd - aq) (it + s)/q, {s, 0, q - 1}]];
   jwr = Flatten[Table[iw (ddd - aq)/q, {q}]]; 
	jt = Flatten[{jtl, jtc, jtr}];
	jw = Flatten[{jwl, jwc, jwr}];		
	nn =  n - If[n <= 2, -1, Floor[n/3]];
      \({{tt, \ ww}, pph, hhi, dd, eer, perr, norm, 
                  ier}\  = \ \ {\ 
                  Transpose[\ 
                    Re[GQW[nn, \ a, \ 
                        b, \ $MachinePrecision\ ]\ ]], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Map[fc, t], Map[fc, jt], \ \ pg, pg, pg, 1, 1}\);
       uu = { .7,  .5,  .2};
       vv = {pg, pg^ .4, pg^ .2} // N;

			  

mmc  =  
                Outer[ kc,  tt,  tt]  .  DiagonalMatrix[ ww];
\(s = \ 
                LinearSolve[\ 
                  IdentityMatrix[\ Length[\ tt]\ ]\  - \ 
                    mmc, 
\ \ \ \ \ \ \ \ Map[\ fc, \ tt\ ]\ ]\);
 mmc1 = 
                Outer[ kc,  jt,  tt ]  .  DiagonalMatrix[ww ];
sh0 =  mmc1  .  s;
 sh  =  sh0  + Map[ fc,  jt ];
 t0 = 
                Flatten[Table[a + (b - a) (it + s)/2, {s, 0, 1}]];
 w0 = Flatten[Table[iw (b - a)/2, {2}]]; 
 qh = 
                Outer[ kc, t0, tt ]  .  DiagonalMatrix[ww ]  .  s + 
                  Map[fc, t0];
 la = 
                sh0 -  Outer[ kc,  jt,  t0 ]  .  
                    DiagonalMatrix[w0 ] . qh;  
 wlist = Map[(wght /. x -> #) &, jt]; If[Min[wlist] < 0, Message[NISolve::wneg]];

     \(If[
                sn < ∞
\ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ eer\  = 
                  If[absrel === 
                      Relative, \ \ \ \((\ 
                          Apply[Plus, wlist\ \ Abs[sh0]^sn\ *jw\ ]\ )\)^\((1/
                            sn)\)/
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \((Apply[Plus, wlist\ \ Abs[sh]^sn\ *jw])\)^\((1/sn)\), \((\ 
                        Apply[Plus, wlist\ \ Abs[sh0]^sn\ *jw\ ]\ )\)^\((1/
                          sn)\)\ ]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ dd = \((Apply[Plus, wlist\ Abs[la]^sn*jw])\)^\((1/
                        sn)\); 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ norm = \((\ 
                      Apply[Plus, wlist\ \ Abs[sh]^sn\ *jw])\)^\((1/
                        sn)\);\ 
\ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ eer\  = 
                  If[absrel === Relative, 
                    Max[wlist\ Abs[sh0]]/
                      Max[wlist\ Abs[sh]], 
Max[
                      wlist\ Abs[
                          sh0]]]; 
\ \ \ \ \ \ \ \ \ dd\  = 
                  Max[wlist\ \ Abs[
                        la]]; 
\ \ \ \ \ \ \ \ \ norm = 
                  Max[wlist\ \ Abs[sh]];
\ \ \ \ \ \ ]\);
		
                                              
 u = { .7,  .5,  .2};
 v = {eer, pg^ .4, pg^ .2} // N;
 e = 1;	
\({{t, \ w}, \ phi, hi, d, err, perr, 
                  norm, ier}\  = \ 
			{\ 
                  Transpose[\ 
                    Re[GQW[\ n, \ a, \ 
                        b, \ $MachinePrecision\ ]\ ]], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ s, sh\ , dd, eer, 10  pg, 
                  norm, 1}\);
		
			   

   Do[  \({{t, \ w}, \ phi, hi, d, err, perr, norm, \ 
                      ier} = \ 
\ \ \ \ \ \ \ NIntegrate2[\ kc, \ fc, \ 
                      x, \ {a, \ b}, \ {t, \ w}, \ phi, ag, pg, wp, 
                      n\ ]; 
\ \ \ \ If[\((\ 
                        ier\  == 0)\) || \((ier == 2)\), \ \ Break[]];\)    
       ,{50} ];
     

 If[ ier == 1, 
                If[absrel === 
                    Relative, 
                    Message[ NISolve::ncvb3,  n , Max[sn0, sn],  
                    PaddedForm[norm, 2], PaddedForm[1.5 perr, 1], 
                    exmp],        
                 Message[ NISolve::ncvb4,  n , Max[sn0, sn],  
                    PaddedForm[norm, 2], PaddedForm[1.5 perr  , 1], 
                    exmp]]];
 If[ ier == 2,  
                If[absrel === 
                    Relative, 
            Message[ NISolve::ncvb3,  n  Round[mxgp/n], Max[sn0, sn], 
                    PaddedForm[norm, 2], PaddedForm[1.5  perr,  1],  
                    exmp],         
            Message[ NISolve::ncvb4,  n Round[mxgp/n], Max[sn0, sn], 
                    PaddedForm[norm, 2], PaddedForm[1.5  perr  , 1], 
                    exmp]]];
     \(If[
                pi === True, 
\ \ \ \ \ Print["L_", 
                  sn0, " norm (with weight ", 
                  exmp, ") of solution -> about ", \ 
                  PaddedForm[norm, 2]]; 
\ \ \ If[
                  absrel === 
                    Relative, 
\ \ \ \ \ Print["Relative L_", sn0, " error (with weight ", 
                    exmp, ") -> about ", \ \ \ \ \ PaddedForm[1.5  perr, 
                      2]]; 
Print["Relative L_", 
                    sn0, " discrepancy (with weight ", 
                    exmp, ") -> about", 
                    PaddedForm[1.3  d, 
                      1]], 
	\ 	\ Print["Absolute L_", 
                    sn0, " error (with weight ", 
                    exmp, ") -> about ", 	\ \ \ \ \ PaddedForm[
                      1.5\ perr\ , 
                      2]]; 
Print["Absolute L_", 
                    sn0, " discrepancy (with weight ", 
                    exmp, ") of solution->about", 
                    PaddedForm[1.3\ \ d, 1]]]
\ \ \ ]\);
		
  

 If[off1, On[LinearSolve::luc]];
     If[off2, On[CompiledFunction::cfn]];
  If[off3, On[General::spell]];
      If[off4, On[NumberForm::sigz]];
            If[off5, On[NIntegrate::inum]];
                 If[off6, On[Set::write]];
                      If[off7, On[NIntegrate::nlim]];

  
		
{{a0,  b0},  kc,  fc, t, wph} = {{a0,  b0},  kc,  
                  fc,  t,  w*phi};
      \(If[\ \((Head[\ kc\ ]\  === \ 
                      CompiledFunction)\)\  && \ \(! MachineNumberQ[\ 
                      x\ ]\), 

                ff[x]\  + \ \((Map[\ \((kk[x, #])\) &, \ t\ ]\  . \ 
                      wph)\), 
\ \ \ \ \ \ fc[
                    x]\  + \ \((Map[\ \((kc[x, #])\) &, \ t\ ]\  . \ 
                      wph)\)]\)
 ];
Not[MatchQ[a, _Complex]  ||  MatchQ[b, _Complex] ];

\(NIntegrate2[
\ kc_, \ fc_, \ 
            x_, \ {\ a_, \ b_}, \ {t_, \ w_}, \ phi_, \ ag_, \ pg_, wp_, \ 
            n_\ ]\  := 
Module[\ 
\ \ \ \ {\ 
\ \  (*\ 
                Local\ variables\ onnected\ with\ Gaussian\ quadratures\ *) 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ t1, 
              kt, ktl, ktc, ktr, w1, kw, kwc, kwl, 
              kwr, 
	
\ \ \  (*\ 
                Connected\ with\ LinearSolve[]*) 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ hh, mc, mm1, 
              ss, 
	
\ \ \  (*\ 
                For\ norms\ of\ solution\ and\ errors\ *) 
\ \ \ \ \ \ \ \ hi1, norm, aer, delta, ts, tst, cm1, cf1, 
              err, 
	
\ \ \  (*For\ automatic\ chosing\ of\ panels*) 
c, c1, c2, c3, r0, r1, r2, pp, pp0, pp1, pp2, 
              pp3, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ff, dff, au, av, cor, 
              q, 
	
\ \  (*\ 
                Remaining\ *) 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ po, ppo, cf2, ap2, aq2
\ \ \ }, 
	

            po = Length[t]/n; 

            c1 = Ceiling[17/n] + po; 
			\  (*\ 
              Solve\ the\ set\ of\ algebraic\ equations\ *) 
	
\ 
            mc = \ Outer[\ kc, \ t, \ t\ ]\  . \ DiagonalMatrix[\ w\ ]; 
\ 
            ss\  = \ LinearSolve[\ 
                IdentityMatrix[\ Length[\ t\ ]\ ]\  - \ mc, 
                Map[\ fc, \ t\ ]\ ]; 

 (*\ 
              Compute\ the\ norm\ of\ the\ relative\ L_sn\ error\ on\ another\ grid\ *) 

\ 
            mm1 = Outer[\ kc, \ jt, \ t\ ]\  . \ DiagonalMatrix[w\ ]; 
\ 
            hi1 = \ mm1\  . \ ss\  + 
                Map[\ fc, \ 
                  jt\ ]; 

\ \ \ \ \ \ \ \ \ If[sn < ∞
\ \ \ \ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ \ \ \ norm = \((\ 
                    Apply[Plus, wlist\ \ \ Abs[hi1]^sn\ *jw])\)^\((1/
                      sn)\)
\ \ \ \ \ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \(norm = 
                  Max[wlist\ Abs[
                        hi1]];\)
\ \ \ ]; 
	


            If[sn < ∞, 
              aer = \((\ 
                    Apply[Plus, wlist\ Abs[hi1\  - \ hi]^sn*jw\ ])\)^\((1/
                      sn)\), aer = 
                Max[wlist\ Abs[hi1\  - \ hi]]]; 	
	
\ \ \ 
\ 
            err\  = If[absrel === Relative, aer\ /norm, aer]; 
\ 
            If[n >= mxgp
, 
perr = 
                err; 
Return[\ {{t, \ w}, \ ss, hi1, d, 
                  err, perr, norm, 1}\ ]\ 
\ ]; 
\ 
            If[err < 1.5  ppr\ 
, 
perr = 
                err; Return[\ {{t, \ w}, \ ss, hi1, d, err, perr, norm, 
                  0}\ ]\ 
\ ]; 		

 (*Compute\ the\ L_sn - discrepancy\ \ and\ a\ prediction\ 
\ \ of\ the\ L_sn\ error\ *) 

            ktc = Flatten[
                Table[a + \((b - a)\) \((it + s)\)/c1, {s, 1, c1 - 2}]]; 

            kwc = Flatten[Table[iw \((b - a)\)/c1, {c1 - 2}]]; 

            ap2 = a + \((b - a)\)/c1; 
q = Ceiling[7/n] + 3; 

            ktl = Flatten[
                Table[a + \((ap2 - a)\) \((it + s)\)/q, {s, 0, q - 1}]]; 

            kwl = Flatten[Table[iw \((ap2 - a)\)/q, {q}]]; \ 

            aq2 = a + \((b - a)\)/c1 \((c1 - 1)\); 

            ktr = Flatten[
                Table[aq2 + \((b - aq2)\) \((it + s)\)/q, {s, 0, q - 1}]]; 

            kwr = Flatten[Table[iw \((b - aq2)\)/q, {q}]]; \ 

            kt = Flatten[{ktl, ktc, ktr}]; 

            kw = Flatten[{kwl, kwc, kwr}]; 		

            hh = Outer[\ kc, kt, t\ ]\  . \ DiagonalMatrix[w\ ]\  . ss\  + \ 
                Map[fc, kt]; 
ts = err/perr; 

            delta = hi1 - \ 
                Outer[\ kc, \ jt, \ kt\ ]\  . \ DiagonalMatrix[\ kw\ ] . 
                  hh - \ \ Map[fc, 
                  jt]; 

\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ If[
              sn < ∞
\ \ \ \ \ \ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(d1 = \((Apply[Plus, 
                        wlist\ Abs[delta]^sn*jw])\)^\((1/
                        sn)\);\)
\ \ \ \ \ \ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(d1 = 
                  Max[wlist\ \ \ Abs[
                        delta]];\)
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 

perr = 
              Min[ .7  err, 
                d1\ err/
                    d]; 

 (*Compute\ the\ correction\ coefficients*) 
				
tst = d1/\((d\ e)\); 
au = Abs[tst - 1]; \ 
            av = Abs[ts - 1]; 
cor = 36\ ArcTan[ .01\ au +  .13  av]/n; 

            If[perr <  .9  
                  ppr\ 
, 
\ \ \ \ \ \ Return[\ {{t, \ 
                    w}, \ ss, hi1, d1, err, perr, norm, 
                  0}\ ]\ ]; 
\ \ \ \ \ \ If[\((po == 
                    np)\) && \((perr >=  .9  ppr)\)\ 
	\ \ , 
	
              Return[\ {{t, \ w}, \ ss, hi1, d1, err, perr, norm, 
                  2}\ ]\ ]; 
\ \ 
 (*Create\ the\ automatic\ choice\ of\ the\ next\ panel' s\ number

                on\ base\ of\ first\ tree\ terms\ of\ the\ error' 
                s\ asymtotic\ expansion*) \ 

\ \ \ \ u = 
              Prepend[Drop[u, \(-1\)], po]; 
v = 
              Prepend[Drop[v, \(-1\)], perr]; 
\ 

            If[po > 2
, \ \ \ \ \ \ 
\ \ \ \ \ pp0 = 
                po*\((perr/\((Abs[1 - u[\([2]\)]/po]\ pg)\))\)^\((1/\((2  
                            n)\))\); 
\ \ \ \ \ cm1 = {{u[\([1]\)]^\((\(-2\) 
                            n)\), u[\([1]\)]^\((\(-2\) n - 
                            2)\)}, \ \ {u[\([2]\)]^\((\(-2\) n)\) - 
                        u[\([1]\)]^\((\(-2\) n)\), 
\ 
                      u[\([2]\)]^\((\(-2\) n - 2)\) - 
                        u[\([1]\)]^\((\(-2\) n - 2)\)}} // 
                  N; 
\ \ \ \ \ cf1 = {v[\([1]\)], 
                  v[\([2]\)]}; \ \ 
\ \ \ \ {c1, c2} = 
                LinearSolve[cm1, cf1]; 
\ \ \ \ ff = 
                c1\ r^\((2  n)\)\  + \ c2\ r^\((2  n + 2)\) - 
                  pg; 
\ \ \ \ dff = D[ff, r]; 
\ \ \ \ r1 = 
                pp0^\((\(-1\))\); 
\ \ \ \ Do[
                r1 = r1 - \((ff /. r -> r1)\)/\((dff /. 
                          r -> r1)\), {8}]; 
\ \ \ \ pp1 = 
                r1^\((\(-1\))\); 
\ \ \ \ cm2 = {{u[\([1]\)]^\((\(-2\) n)\), 
                      u[\([1]\)]^\((\(-2\) n - 2)\), 
                      u[\([1]\)]^\((\(-2\) n - 
                            4)\)}, 
\ \ \ \ \ \ \ \ \ \ {u[\([2]\)]^\((\(-2\) n)\) - u[\([1]\)]^\((\(-2\) n)\), 
                      u[\([2]\)]^\((\(-2\) n - 2)\) - 
                        u[\([1]\)]^
		\((\(-2\) n - 2)\), 
                      u[\([2]\)]^\((\(-2\) n - 4)\) - 
                        u[\([1]\)]^\((\(-2\) n - 
                              4)\)}, 
\ \ \ \ \ \ \ \ \ \ {u[\([3]\)]^\((\(-2\) n)\) - u[\([2]\)]^\((\(-2\) n)\), 
                      u[\([3]\)]^\((\(-2\) n - 
                              2)\) - 
\ \ \ \ \ \ \ \ \ \ \ u[\([2]\)]^\((\(-2\) n)\), 
                      u[\([3]\)]^\((\(-2\) n - 4)\) - 
                        u[\([2]\)]^\((\(-2\) n)\)}} // N; 

              cf2 = {v[\([1]\)], v[\([2]\)], v[\([3]\)]}; 
{c1, c2, c3} = 
                LinearSolve[cm2, cf2]; 

              ff = c1\ r^\((2  n)\)\  + \ c2\ r^\((2  n + 2)\) + 
                  c3\ r^\((2  n + 4)\) - pg; 
dff = D[ff, r]; 

              r0 = \((po*\((perr/\((Abs[
                                  1 - u[\([2]\)]/po]\ pg)\))\)^\((1/\((2  
                                n)\))\))\)^\((\(-1\))\); 

              r2 = pp1^\((\(-1\))\); 

              Do[r2 = r2 - \((ff /. r -> r2)\)/\((dff /. r -> r2)\), {6}]; 

              pp2 = r2^\((\(-1\))\);
\ \ ]; 

            pp3 = Which[po < 4, pp0, po < 5, pp1, po < np + 1, pp2]; 

            pp = If[po > 
                  2
\ \ \ \ \ \ \ , 
	\ \ \ \ \ \ Round[
                  pp3 \((1 + cor/17)\) + cor/5.  + 1. /n], 
                po + 1]; 

		 (*\ 
              Prediction\ of\ the\ next\ panel' s\ number\ *) 
		\ 

            pp = Max[po + 1, pp]; 
\ 
		 (*\ 
              Scale\ the\ \ new\ quadrature\ for\ the\ next\ step*) 


            If[po < 2 + \((25 + pp/10)\)/
                    n
, 
\ \ \ \ \ c = \ 
                Min[pp, 
                  Ceiling[\((1 + 2. /n)\)\ po]]
\ \ \ \ \ \ , 
\ \ \ \ \ c := 
                Min[pp, 
                  np]]; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ If[\((c <= po)\), c := Min[po + 1, np]]; \ 
\ \ \ \ \ If[c > 3, 
              q := 2/\((np - po + 1)\), q := 0]; \ 
\ \ \ \ \ ppr = 
              pg\ l^q; 
\ \ \ \ \ If[
              perr <  .9\ ppr
\ \ \ \ , \ \ \ 
Return[\ {{t, \ w}, \ ss, hi1, d1, err, perr, norm, 
                  0}]]; \ \ \ 
\ \ \ \ \ \ \ If[c >  .9\ np, 
              c := np]; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ t1 = 
              Flatten[Table[
                  a + \((b - a)\) \((it + s)\)/c, {s, 0, 
                    c - 1}]]; 
\ \ \ \ \ \ \ w1 = 
              Flatten[Table[
                  iw \((b - a)\)/c, {c}]]; 
\ \ \ \ \ \ \ e = \((po/
                    c)\)^\((2  n)\); 
\ 
			 (*\ 
              Return\ the\ new\ grid\ and\ solution\ on\ it\ \ *) 
\ \ 
\ \ \ Return[\ {{t1, \ w1}, \ ss, hi1, d1, err, perr, norm, 1}\ ]; 
\ \ \ Return[
              u, \ v, \ ppr, e, asa, asb];
\ \ \ ];\)



\(GQW[n_, \ a_, \ b_, \ 
            prec_]\  := 

 (*
\ \ \ Reproduced\ without\ modifications\ from\ \(package

              \ \ \ \ \ \ \(\(:\)\(Name : \ \(NumericalMath`GaussianQuadrature`

                  \ \ \ \ \ \ \(\(:\)\(Title : \ 
                      Coefficients\ for\ Gaussian\ \(Quadrature

                        \ \ \ \ \ \ \(\(:\)\(Author : \ 
                            Jerry\ B . \ \(Keiper

                                \ \ \ \ \ \ \(\(:\)\(Package\ \(Version : \ \(Mathematica\ 2.0

                                      \ \ \ \ \ \ \(\(:\)\(Copyright : \ 
                                        Copyright\ 1990\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\), 
\ \ \ \ \ \ Wolfram\ Research, \(\(Inc\)\(.\)\)
\ *) 
\ \ \ \ \ \ 

          Module[{i, \ w, \ x, \ m, \ rhs, \ legn, \ 
              t}, 
\ \ \ \ \ \ \ \ \ legn\  = \ 
              Numerator[LegendreP[n, t]]; 
	
            If[\ n\  == \ 
                1
\ \ \ \ \ \ \ \ , 
\ \ \ \ \ \ \ \ Return[N[{{0, \ 2}}, \ prec]]]; 
	
            x\  = \ Apply[List, \ 
                Map[#[\([2]\)] &, \ NRoots[legn\  == \ 0, \ t, \ prec]]]; 
	
            m\  = \ Table[x^i, \ {i, n - 1}]; 
            PrependTo[m, \ Table[1, \ {n}]]; 
	
            rhs\  = \ 
              Table[\((1 + Sign[\((\(-1\))\)^i])\)/\((i + 1)\), {i, 0, 
                  n - 1}]; 
	
            w\  = \ Expand[LinearSolve[m, \ rhs] \((b - a)\)/2]; 
	
            x\  = \ Expand[\((x\ \((b - a)\)\  + \ \((b + a)\))\)/2]; 
	
            Transpose[{x, \ w}]
		];\)



\(NISolveVolterra[\ {k_, var_, aa_, bb_, f_}, \ {var0_, \ a0_, \ 
              b0_}, \ \ opt___\ ] := 
	
          Module[{fi, ww, pg, ll, llpos, aanew, bbnew, a0new, b0new, knew, 
              fnew, sol, 
              ss}, 
{fi, ww, 
                pg} = \({SymbolicSolution, Weight, PrecisionGoal} /. {opt}\) /. 
                Options[NISolveVolterra]; 
acont = 
              a0; 
bcont = b0; 
acont0 = 
              aa; 
bcont0 = 
              bb; 
fintopt = 
              FilterOptions[FunctionInterpolation, 
                opt]; 
If[\((fi === True)\) || \((fi === 
                    False)\), , Message[NISolve::symsol]; 
              fi = False]; 
exmp = 
              If[ww === Automatic, 1, ww]; 
If[
              ww === Automatic, ww = 1]; 
If[\ 
              pg\  === \ Automatic, \ pg\  = \ 5, 
\ 
              If[Not[IntegerQ[pg]] || Not[Positive[pg]], 
                Message[NISolve::prg]; 
                pg = $MachinePrecision - 
                    11]]; 

ll = 
              Flatten[{aa, bb, a0, b0}]; 
llpos = 
              Flatten[{Position[ll, ∞], 
                  Position[ll, \(-∞\)]}]; 
	
            If[Length[llpos] > 
                0, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ aanew = 
                ArcTan[aa /. 
                    var0 -> Tan[var0]]; 
	\ \ \ 		\ \ bbnew = 
                ArcTan[bb /. 
                    var0 -> 
                      Tan[var0]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ a0new = 
                ArcTan[a0]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ b0new = ArcTan[b0]; 
		\ \ \ 	\ \ knew = 
                Table[\((k /. {var0 -> Tan[var0], 
                            var[\([i]\)] -> Tan[var[\([i]\)]]})\)[\([i]\)]/
                    Cos[var[\([i]\)]]^2, \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {i, 1, 
                    Length[k]}]; 
	\ \ 	\ \ fnew = \((f /. 
                    var0 -> Tan[
                        var0])\); 
\ \ \ \ \ \ \ \ \ \ \ \ \ If[ww === Automatic, 
	\ \ \ \ sol = 
                  NISolveVolterraLimit[{knew, var, aanew, bbnew, 
                      fnew}, {var0, \ a0new, \ b0new}, opt, \ 
                    Weight -> 
                      1/Cos[var0]^2]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ If[fi, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ss[
                      z_] = sol[
                      ArcTan[z]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Format[ss] := 
                    SequenceForm["SymbolicSolution[", {a0, 
                        b0}, ",<>]"]; \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                    ss]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ , \ \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ sol = FunctionInterpolation[Evaluate[sol[x]], {x, ArcTan[a0], ArcTan[b0]}, 
                      Evaluate[
                        fintopt], 
\ \ \ \ \ \ \ \ \ \ \ InterpolationPrecision -> 
                        Min[pg + 
                            1, $MachinePrecision], 
\ \ \ \ \ \ \ \ \ \ \ \ InterpolationOrder -> 7, 
                      PrecisionGoal -> Min[pg + 1, $MachinePrecision], 
                      AccuracyGoal -> 
                        Min[$MachinePrecision, 
                          pg + 1], \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ InterpolationPoints -> 200, 
                      MaxRecursion -> 
                        50]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ss[z_] = sol[ArcTan[
                        z]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Format[ss] := SequenceForm["InterpolatingFunction[", {{a0, 
                          b0}}, ",<>]"]; \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                    ss]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ], 
sol = 
                  NISolveVolterraLimit[\ {knew, \ var, aanew, bbnew, 
                      fnew}, {var0, \ a0new, \ b0new}, \ 
                    Weight -> \((ww /. var0 -> Tan[var0])\)/Cos[var0]^2, 
                    opt]; 
\ \ \ \ \ \ \ \ \ \ \ \ If[
                  fi, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ss[z_] = 
                    sol[ArcTan[
                        z]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Format[ss] := 
                    SequenceForm["SymbolicSolution[", {a0, 
                        b0}, ",<>]"]; \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                    ss]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ , \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ sol = 
                    FunctionInterpolation[
                      Evaluate[sol[x]], {x, a0new, b0new}, \ 
                      Evaluate[
                        fintopt], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ InterpolationPrecision -> 
                        Min[pg + 
                            1, $MachinePrecision], 
\ \ \ \ \ \ \ \ \ \ \ \ InterpolationOrder -> 7, 
                      PrecisionGoal -> Min[pg + 1, $MachinePrecision], 
                      AccuracyGoal -> 
                        Min[$MachinePrecision, 
                          pg + 1], \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ InterpolationPoints -> 200, 
                      MaxRecursion -> 
                        50]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ss[z_] = sol[ArcTan[
                        z]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Format[ss] := SequenceForm["InterpolatingFunction[", {{a0, 
                          b0}}, ",<>]"]; \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                    ss]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]
\ \ \ \ \ \ \ \ \ \ \ ]
, 
sol = 
                NISolveVolterraLimit[\ {k, var, aa, bb, f}, \ {var0, \ a0, \ 
                    b0}, \ opt\ ]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ If[fi, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ss[
                    z_] = sol[
                    z]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Format[ss] := SequenceForm["SymbolicSolution[", {a0, 
                      b0}, ",<>]"]; \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                  ss]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ , \ \ \ \ \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ sol = 
                  FunctionInterpolation[Evaluate[sol[x]], {x, a0, b0}, 
                    Evaluate[
                      fintopt], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ InterpolationPrecision -> 
                      Min[pg + 
                          1, $MachinePrecision], 
\ \ \ \ \ \ \ \ \ \ \ \ InterpolationOrder -> 7, 
                    PrecisionGoal -> Min[pg + 1, $MachinePrecision], 
                    AccuracyGoal -> 
                      Min[$MachinePrecision, 
                        pg + 1], \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ InterpolationPoints -> 200, 
                    MaxRecursion -> 
                      50]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ss[z_] = 
                  sol[z]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Format[ss] := SequenceForm["InterpolatingFunction[", {{a0, 
                        b0}}, ",<>]"]; \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Return[
                  ss]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]

]
		];\)


NISolveVolterraLimit[{kg_List, var_List, aa1_List, bb1_List, 
            f1_}, {var0_, ag_, bg_}, opt___]:=Module[{wght,pg,mxgp,pi,sn,absrel,var1,mm,kold,k,knew,f1f,gf,g,aa1func,bb1func,cf,c,cs,c1,c1f,icf,ic,aaf,bbf,aac,bbc,a1f,f,a1,b1,flst,coef,solonlSlow, solonlFast,solonllwithlSlow,solonllwithlFast,nnk,lnnk,l1,l2,nnc,rrnnc,tc,tcunf,iop,wlist,sh,nn,t,solt,w,err,step,nndo,tdo,soltdo,wdo,errdo,norm,cflist,c1flist,bnew,b,kL},

  

 {wght, pg, mxgp, pi, sn, 
                  absrel} = {Weight, PrecisionGoal, MaxGridPoints, 
                      PrecisionInformation, SolutionNorm, Error} /. {opt} /. 
                  Options[NISolveVolterra];


If[wght === Automatic, wght = 1];
If[(pi === True) || (pi === 
                      False), , Message[NISolve::precinf]; pi = False];	
If[ pg  ===  Automatic,  
                pg  =  5, 
 
                If[Not[IntegerQ[pg]] || Not[Positive[pg]], 
                  Message[NISolve::prg]; pg = 5]];
pg = (10^(-pg));
If[mxgp === Automatic, mxgp = 200];
If[
                Not[IntegerQ[mxgp]] || Not[Positive[mxgp]], 
                Message[NISolve::mxgpmes]; mxgp = 200];

If[(absrel === Relative) || (absrel === 
                      Absolute), , Message[NISolve::errmes]; 
                absrel = Relative];
   If[(sn === ∞) || (sn === 
                      Automatic), , 
                If[Not[IntegerQ[sn]] || Not[Positive[sn]], 
                  Message[NISolve::nrmes]; sn = 2]];
If[sn === Automatic, sn = 2];

var1 = Prepend[var, var0];
mm = Length[kg];
kold = (bg - ag) kg /. 
                  Flatten[
                    Transpose[{var1, var1}] /. {i_Symbol, 
                          i_Symbol} -> {i -> ag + (bg - ag) i}];
       k = 
                Table[Compile[Evaluate[{var0, {var[[i]], _Real, 1}}], 
                    Evaluate[kold[[i]]]], {i, 1, mm}];
knew = 
                Table[Function[Evaluate[{var0, var[[i]]}], 
                    Evaluate[kold[[i]]]], {i, 1, mm}];

f1f = Function[{var0}, Evaluate[f1]];
gf[x_] = f1f[ag + (bg - ag) x];
SetAttributes[gf, Listable];


aa1func = 
                Function[{var0}, Evaluate[aa1]];
bb1func = 
                Function[{var0}, Evaluate[bb1]];

cf = 
                Function[{x}, Evaluate[Sin[π  x/2]^2]];
c = Compile[{x}, Evaluate[cf[x]]];
cs = 
                Compile[{{x, _Real, 1}}, Evaluate[cf[x]]];
c1 = 
                Compile[{{x, _Real, 1}}, Evaluate[ Sin[π x]]]; 
c1f = 
                Function[{x}, Evaluate[ Sin[π x]]]; 
icf = 
                Function[{x}, ArcCos[1 - 2 x]/π];
ic = Compile[{x}, Evaluate[icf[x]]];
aaf = 
                Function[{x}, 
                  Evaluate[(1/(bg - 
                            ag)) (aa1func[(ag + (bg - ag) x)] - 
                          ag)]];
bbf = 
                Function[{x}, 
                  Evaluate[(1/(bg - 
                            ag)) (bb1func[(ag + (bg - ag) x)] - 
                          ag)]];
aac = Compile[{x}, Evaluate[aaf[x]]];
bbc = Compile[{x}, Evaluate[bbf[x]]];
a1f = 
                Function[{x}, Evaluate[icf[aaf[cf[x]]]]];
b1f = 
                Function[{x}, Evaluate[icf[bbf[cf[x]]]]];
a1 = 
                Compile[{x}, Evaluate[icf[aaf[cf[x]]]]];
b1 = 
                Compile[{x}, Evaluate[icf[bbf[cf[x]]]]];


flst[lf_]:=Block[{fy, m},
m = Length[lf];
fy = RotateLeft[lf, Floor[m/2]];
RotateRight[Fourier[fy]/(Sqrt[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, 
                    PeriodicInterpolation -> True]
];

 solonlSlow[l_List]:=Module[{nn, clist, bnew, bb, b, kL, aux, mat, ff},
nn = Length[l];
clist = cs[l];
bnew[x_, 
                            z_] = coef[nn][z - x] + coef[nn][z + x];
bb[
                            x_] = ((1  )/((2  nn + 1. ))) 
                            Map[bnew[a1f[x], #] - bnew[b1f[x], #] &, l];
b = Transpose[Map[bb[#] &, l]];
kL = 
                          Transpose[
                            Map[Table[(k[[i]][#, cs[l]]*
                                      c1[l]), {i, 1, mm}] &, clist], {2, 3, 
                              1}];
aux = Apply[Plus, b   kL, {2}];
 mat = IdentityMatrix[Length[l]] - aux;
ff = gf[clist];
LinearSolve[Transpose[mat], ff]
];


 solonlFast[l_List]:=Module[{nn, clist, an, bn, pop, rr, b, kL, aux, mat, 
                        ff},
nn = Length[l];
clist = cs[l];
an = Map[a1[#] &, l];
bn = Map[b1[#] &, l];
pop = 
                          Compile[{n, {x, _Real, 1}}, 
                            Evaluate[(1/n) (Cos[π n an] - 
                                    Cos[π n bn]) Sin[π n x]]];
rr = Range[1, nn];
b = (2. /((2  nn + 1. ))) 
                            Transpose[
                              Apply[Plus, 
                                Transpose[Map[pop[#, l] &, rr]], {1}], {2, 3, 
                                1}];
kL = 
                          Transpose[
                            Map[Table[(k[[i]][#, cs[l]]*
                                      c1[l]), {i, 1, mm}] &, clist], {2, 3, 
                              1}];
aux = Apply[Plus, b   kL, {2}];
 mat = IdentityMatrix[Length[l]] - aux;
ff = gf[clist];
LinearSolve[Transpose[mat], ff]
];

 
solonllwithlSlow[l_List, ll_List, yy_List]:=Module[{clist, an, bn, b, kL, ff, nn1, rr, pop, bnew, 
                        bb},
nn1 = Length[l];
bnew[x_, 
                            z_] = coef[nn1][z - x] + coef[nn1][
                              z + x];
bb[
                            x_] = ((1  )/((2  nn1 + 1. ))) 
                            Map[bnew[a1f[icf[x]], #] - 
                                  bnew[b1f[icf[x]], #] &, l];  
b = Transpose[Map[bb[#] &, ll]];
kL = 
                          Transpose[
                            Map[Table[(k[[i]][#, cs[l]]*
                                      c1[l]), {i, 1, mm}] &, ll], {2, 3, 
                              1}];
ff = gf[ll];
ff + 
                          Apply[Plus, Transpose[yy . (b   kL)]]
];


solonllwithlFast[l_List, ll_List, yy_List]:=Module[{clist, an, bn, b, kL, ff, nn1, rr, pop, bnew, 
                        bb},
nn1 = Length[l];
clist = Map[ic[#] &, ll];
an = Map[a1[#] &, clist];
bn = Map[b1[#] &, clist];
pop = 
                          Compile[{n, {x, _Real, 1}}, 
                            Evaluate[(1/n) (Cos[π n an] - 
                                    Cos[π n bn]) Sin[π n x]]];
 rr = Range[1, nn1];
 b = (2. /((2  nn1 + 1. ))) 
                            Transpose[
                              Apply[Plus, 
                                Transpose[Map[pop[#, l] &, rr]], {1}], {2, 3, 
                                1}];
kL = 
                          Transpose[
                            Map[Table[(k[[i]][#, cs[l]]*
                                      c1[l]), {i, 1, mm}] &, ll], {2, 3, 
                              1}];
ff = gf[ll];
ff + 
                          Apply[Plus, Transpose[yy . (b   kL)]]
];

nnk = 
                Round[mxgp (bg - ag) // N];
lnnk = Table[k/nnk, {k, 0, nnk}];
l1 = 
                Sort[Flatten[{acont0 /. (var0 -> 
                                ag + (bg - ag) lnnk)}] // N] // First;
l2 = 
                Sort[Flatten[{bcont0 /. (var0 -> 
                                ag + (bg - ag) lnnk)}] // N] // Last;
Which[(acont == -∞) && (bcont == ∞), , (acont == -∞) && ( bcont < l2), 
                Message[NISolve::limitincons, acont, bcont]; 
                Break[$Failed], (bcont == ∞) && (l1 <  
                      acont), Message[NISolve::limitincons, acont, bcont]; 
                Break[$Failed], (l2 > bcont) || (l1 <  acont), 
                Message[NISolve::limitincons, acont, bcont]; 
                Break[$Failed]];

nnc = Round[ mxgp/4];
rrnnc = Range[1, nnc];
tc = 
                Table[Sin[π/2   2.  s/(2 nnc + 1)]^2, {s, 1, 
                    nnc}];
tcunf = 
                Table[ 2.  s/(2 nnc + 1), {s, 1, nnc}];
\(iop = \ 
                Sin[π\ \ tcunf]\ \ Apply[Plus, 
                    Map[\(\(1.  - \((\(-1\))\)\^#\)\/#\) 
                          Sin[π\ #\ \ tcunf] &, rrnnc]]\);
wlist = 
                Table[(wght /. var0 -> tc[[i]]), {i, 1, 
                    Length[tc]}];
If[Min[wlist] < 0, 
                Message[NISolve::wneg]];

sh = 120;

nn = 1;
t = 
                Table[2.   s/(2 nn + 1), {s, 1, nn}];
solt = 
                If[nn ≤ sh, solonlFast[t], solonlSlow[t]];
w = 
                If[nn ≤ sh, solonllwithlFast[t, tc, solt], 
                  solonllwithlSlow[t, tc, solt]];
err = pg (10^11);

step = 1;

  Do[
 
                Which[ 
           step == 
                    1, 
                               nndo =  
                    Which[   (10^10) <= 
                        err/pg,                    10 + 
                        Round[1.5  nn], 
                                                       (10^8) <= err/pg < (10^10),           10 +  
                        Round[1.3  nn], 
                                                        (10^6) <= err/pg < (10^8),           12 +  
                        Round[1.2  nn], 
                                                        (10^4) <= err/pg < (10^6),           13 + 
                        Round[1.1  nn], 
                                                        (10^2) <= err/pg < (10^4),           12 + 
                        Round[1.05  nn], 
                                                         err/pg < 
                        (10^2),                        nn + 10, 
                                                         err/pg < 
                        10 ,                         nn + 5 ], 
           step == 
                    2, 
                           nndo =  
                    Which[   (10^10) <= 
                        err/pg,                  3 + 
                        Round[1.1  nn], 
                                                    (10^8) <= err/pg < (10^10),         3 + 
                        Round[1.05  nn], 
                                                    (10^6) <= err/pg < (10^8),           3 + 
                        Round[1.02  nn], 
                                                    (10^4) <= err/pg < (10^6),           3 + 
                        Round[1.01  nn], 
                                                    (10^2) <= err/pg < (10^4),           nn + 
                        5, 
                                                     err/
                          pg < (10^2),                        nn + 
                        4, 
                                                     err/
                          pg < 10 ,                         nn + 3 ]
]; 
If[
                  nndo < mxgp - 0.4 mxgp, step = 1, 
                  step = 2]; 
tdo = 
                  Table[2.  s/(2  nndo + 1. ), {s, 1, 
                      nndo}]; 
soltdo = 
                  If[nndo ≤ sh, solonlFast[tdo], 
                    solonlSlow[tdo]]; 
wdo = 
                  If[nndo ≤ sh, solonllwithlFast[tdo, tc, soltdo], 
                    solonllwithlSlow[tdo, tc, soltdo]]; 
If[
                  sn < ∞, 
             errdo = (((Apply[Plus, (2/(2 nnc + 1)) wlist  Abs[wdo - w]^sn  iop]))^(1/sn)); 
             norm = (((bg - ag))^(1/
                            sn)) (((Apply[Plus, (2/(2 nnc + 1)) wlist  Abs[wdo]^sn  iop]))^(1/sn)); 
             If[absrel === Relative, 
                    errdo = 
                      errdo/norm
                 ], 
             errdo = 
                    Max[wlist Abs[
                          wdo - w]]; 
             norm = Max[wlist  Abs[
                          wdo]]; 
            If[absrel === Relative,   errdo = 
                      errdo/norm];
     ]; 
If[errdo ≤ 3 pg, s2 = nn; s3 = t; s4 = solt; 
                  err = errdo; Return[$Failed]]; 
If[
                  nndo ≥  mxgp, s2 = nn; s3 = t; s4 = solt; 
                  err = errdo; 
     If[
                    absrel === Relative,  
                    Message[ NISolve::ncvb,  s2 , sn,  
                      PaddedForm[norm, 2], PaddedForm[ err, 1], 
                      exmp], 
                                   Message[ NISolve::ncvb2, s2 , sn,  
                      PaddedForm[norm, 2], PaddedForm[ err, 1], 
                      exmp]]; 
      Return[$Failed]
]; 
 
                nn = nndo; 
t = 
                  tdo; 
solt = 
                  soltdo; 
w = 
                  wdo; 
err = errdo,   {(10^8)}];


 \(If[
                pi === True, 
Print["L_", 
                  sn, " norm  (with weight ", 
                  exmp, ") of solution -> about ", \ 
                  PaddedForm[norm, 3]]; 
\ 
                If[
absrel === 
                    Relative, \ \ \ \ \ 
Print["Relative L_", sn, " error (with weight ", exmp, ") -> about ", 
                    PaddedForm[\ err, 
                      3]]
\ \ \ \ \ \ \ \ , 
Print["Absolute L_", sn, " error (with weight ", 
                    exmp, ")-> about ", \ 
                    PaddedForm[\ err, 
                      3]]
\ \ \ \ \ ]
\ \ \ \ \ ]\);

cflist = cs[s3];
c1flist = Sin[π  s3] // N;
bnew[x_, 
                  z_] = coef[s2][z - x] + coef[s2][z + x];
b[
                  x_] = ((1  )/((2  s2 + 1. ))) 
                  Map[bnew[a1f[icf[x]], #] - bnew[b1f[icf[x]], #] &, s3];
kL[x0_]=Transpose[Table[(knew[[i]][x0, cflist]*
                          c1flist), {i, 1, mm}]];
Function[{var0},Evaluate[f1f[var0]+Apply[Plus, 
                          s4 . (b[(var0 - ag)/(bg - 
                                        ag)]  kL[(var0 - ag)/(bg - 
                                        ag)])]]]

];
Options[NISolveFredholm] = 
        Drop[Options[
            NISolveFredholm], {3}];
Options[
          NISolveVolterra] = 
        Drop[Options[NISolveVolterra], {2}];
Options[NISolve]={{"FOR FREDHOLM TYPE EQUATIONS ->", Options[NISolveFredholm]},{"FOR VOLTERRA TYPE EQUATIONS ->", Options[NISolveVolterra]}}; 
End[];
EndPackage[]

Introduction to NISolve[ ]

The package solves numerically an integral equation of the form

y (x) = 
          f (x) + (∫_a^b) 
                k (x, t) y (t) \[DifferentialD]t, a≤x≤b                 (1)

with some constants a and b , -∞≤a, b≤∞, when unique solvability is provided in the prescribed weighted (Subscript[L, p]) - spaces (i.e. corresponding  integral operator is invertable).  The algorithm is the best when both the kernel k(x,t) and the free term f(x) are  rather smooth and have not sharp changing derivatives (i.e. have no high oscillations, peaks etc.).  In other cases the algorithm works slowly and sometimes cannot provide the prescribed precision.  An efficient solution of equation (1) is provided also when the kernel  k(x,t) is not smooth on some lines in  the domain (a,b)x(a,b) of (x,t)-plane but we can rewrite (1) in the form

y (x) = f (x) + (∫_((Subscript[a, 1]))^((Subscript[b, 1])))(Subscript[k, 1]) (x, 
                    t) y (t) \[DifferentialD]t +  ...  + (∫_((Subscript[a, n]))^((Subscript[b, n])))(Subscript[k, n]) (x, 
                t) y (t) \[DifferentialD]t        (2)

with rather smooth free term f(x),kernels (Subscript[k, p]) (x, t)and (Subscript[a, p]) = (Subscript[a, p]) (x),(Subscript[b, p]) = (Subscript[b, p]) (x), p=1,...,n. For example, it provides the solution of Volterra type equations and boundary problems for differential-functional equations (see below Examples).  
   Taking into account the described circumstances the package is supported on two completely different algorithms. One of them provides solution of (1) using Gauss quadratures. This solver was developed for WRI in August 1997 and now is essentially improved. The second algorithm was developed during 1998-2002 and provides solution of (2) using the generalized Euler-Maclaurin formula introduced in an article of the author of this package. 
   The output of NISolve[] may be InterpolatingFunction (default) or SymbolicSolution object. The latter (called earlier as a NaturalInterpolatingFunction) has the form of a functional sum and (as a rule) is more fast and precise but sometimes rather complicated.In the first 
algorithm precision control is founded both on recursive error and discrepancy in prescribed spaces. The discrepancy for (1) mean the norm of the function

dis (x)=(Overscript[y, n]) (x)-f (x)-(∫_a^b)k(x, t) (Overscript[y, n])(t)\[DifferentialD]t                 (3)

when (Overscript[y, n]) (x) is approximate solution in  n - 
          th step of the recursion . It may be called using  the 
option  SolutionInformation->True.

NISolve[eqns,y,{x,(Subscript[x, min]),(Subscript[x, max])}]       \(Solve\ numericaly\ for\ the\ function\ y, 
\(\(with\)\(\ \)\(the\)\(\ \)\(independent\)\(\ \)\(variable\)\(\ \)\(x\)\(\ \)\(in\)\(\ \)\)\)
the range xmin to xmax

Numerical solution of one - 
        dimensional linear integral equations of the second kind

In addition to all FunctionInterpolation options the following options can be given

Options[FunctionInterpolation]

Option name                           Option   value                                                                                                             

SymbolicSolution
\(\(Compiled\^*\)\)
\(\(GaussPoints\^*\)\)
MaxGridPoints
PrecisionGoal
PrecisionInformation
SolutionNorm
Weight
Error      False
True
Automatic
Automatic
Automatic
False
Automatic
Automatic
Relative              form of output
whether to compile the original equations
number of Gauss points on each panel
total number of grid points
digits of precision sought
information about error and norm of solution
number of p in (Subscript[L, p]) space (p = ∞ is a uniform norm)
weight in  (Subscript[L, p]) space
absolute or relative error to use

*This option is available only for Fredholm  type  equations when  (Subscript[a, 1]) =.  ..  = (Subscript[a, n]) = const,  (Subscript[b, 1]) =.  ..  = (Subscript[b, n]) = const.

1. Easy Examples

1.1

Consider the following equation

y (x) - (∫_1^3) 
                  (sin^2) (x t) 
              y (t) \[DifferentialD]t = ((sin (6 x) - sin (2 x))/(4 x))

with exact solution y (x) = 1.

We construct a  high-precision solution with relative (Subscript[L, 2]) (1, 3)-error  and output give as  InterpolatingFunction object. Besides we demand from NISolve to show information about the error, discrepancy and norm of the solution.

sol = 
      NISolve[y[
                x] - (∫_1^3) Sin[x t]^2 
                  y[t] \[DifferentialD]t == ((Sin[6 x] - Sin[2 x])/(4 x)), 
          y, {x, 1, 3}, 
PrecisionGoal -> 
            14, 
PrecisionInformation -> 
            True
] // Timing

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

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

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

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

Here is the actual (Subscript[L, 2]) discrepancy of the solution

NIntegrate[Abs[sol[[2]][x] - NIntegrate[Sin[x t]^2 sol[[2]][t], {t, 1, 3}] - ((Sin[6 x] - Sin[2 x])/(4 x))]^2, {x, 1, 3}]^(1/2)

1.2

Consider an equation with infinite limit of integration

y (x) - (∫_2^∞) 
            y (t) (E^(-x  t)) \[DifferentialD]t = 
      (e^(-x)) - ((e^(-2 (1 + x)))/(1 + x))

with exact solution y (x) = (e^(-x)).

We construct a high-precision solution with (Subscript[L, 3]) (2, ∞)-error and (Subscript[L, 3]) (2, ∞)-norm of the solution

sol = 
      NISolve[y[
                x] - (∫_2^∞) y[
                    t] (E^(-x  t)) \[DifferentialD]t == 
            (E^(-x)) - ((E^(-2 (1 + x)))/(1 + x)), 
          y, {x, 2, ∞}, 
PrecisionGoal -> 
            14, 
SolutionNorm -> 
            3, 
Error -> 
            Absolute, 
PrecisionInformation -> True] // 
        Timing

Here is  the actual  absolute (Subscript[L, 3])-error

NIntegrate[Evaluate[ Abs[sol[[2]][x] - (E^(-x))]^3], {x, 2, ∞}]^(1/3)

Here is the actual  (Subscript[L, 3]) - norm of the solution

NIntegrate[Evaluate[ Abs[(E^(-x))]^3], {x, 2, ∞}]^(1/3)

Here is the actual  (Subscript[L, 3]) - discrepancy

NIntegrate[Abs[sol[[2]][x] - NIntegrate[(E^(-x t)) sol[[2]][t], {t, 2, 200}] - ((E^(-x)) - ((E^(-2 (1 + x)))/(1 + x)))]^3, {x, 2, 200}]^(1/3)

1.3

This tries to generate a high-precision solution to a complex-valued equation. The exact solution is unknown.

sol = 
      NISolve[(E^x) + I x^3 - (E^(-2 x)) y[x] + 
              Sin[3  x - 
                    2  I] (∫_(-2)^1.5) y[t] Sin[
                      t - I x (E^t)] \[DifferentialD]t == 0, 
          y, {x, -2, 1.5}, 
PrecisionGoal -> 
            12, 
PrecisionInformation -> True] // 
        Timing

This uses more GaussPoints and reaches to prescribed precision faster

sol = 
      NISolve[(E^x) + I x^3 - (E^(-2 x)) y[x] + 
              Sin[3  x - 
                    2  I] (∫_(-2)^1.5) y[t] Sin[
                      t - I x (E^t)] \[DifferentialD]t == 0, 
          y, {x, -2, 1.5}, 
PrecisionGoal -> 
            12, 
GaussPoints -> 
            15, 
PrecisionInformation -> True] // 
        Timing

Here is graps of real and imanginery parts of the solution.

Plot[{Re[sol[[2]][x]], Im[sol[[2]][x]]}, {x, -2, 
        1.5}, PlotStyle -> {RGBColor[0.996109,  0,  0], 
          RGBColor[0,  0,  0.996109]}]

1.4

Consider the follwing Volterra type integral equation

y (x) +  (∫_0^x) (Subscript[J, 0]) (x - t) 
            y (t) \[DifferentialD]t = 
      sin (x)  + (Subscript[J, 0]) (x)

with exact solution y(x)=(Subscript[J, 0])(x).

We construct a high-precision solution with absolute uniform error

sol = 
      NISolve[ y[
                x] +  (∫_0^x) BesselJ[0, x - t] 
                  y[t] \[DifferentialD]t == Sin[x]  + BesselJ[0, x], 
          y, {x, 0, 1}, 
PrecisionGoal -> 
            14, 
Error -> 
            Absolute, 
SolutionNorm -> ∞
] // Timing

Here is the actual absolute error

Plot[Evaluate[Abs[sol[[2]][x] - BesselJ[0, x]]], {x, 0, 1}, 
      PlotRange -> All]

1.5

Consider the following equation

y[x] - 
        3 (∫_1^(x + 1))(1/(((t + x))^2)) 
              y[t] \[DifferentialD]t == (((E^(-1 - x)) (3 (1 + x) - 3 (E^x) (1 + 2 x) + E (1 + 3 x + 2 (x^2))))/((1 + x) (1 + 2 x))) + 3 (E^x) (ExpIntegralEi[-1 - 2 x] - 
              ExpIntegralEi[-1 - x])

with  exact solution y (x) = (e^(-x)).We increase the InterpolationPoints for more exact  output (default=200).

sol = 
      NISolve[y[x] - 
              3 (∫_1^(x + 1))(1/(((t + x))^2)) 
                    y[t] \[DifferentialD]t == (((E^(-1 - x)) (3 (1 + x) - 3 (E^x) (1 + 2 x) + E (1 + 3 x + 2 (x^2))))/((1 + x) (1 + 2 x))) + 
              3 (E^x) (ExpIntegralEi[-1 - 2 x] - 
                    ExpIntegralEi[-1 - x]), 
y, {x, 
            1, 3}, 
PrecisionGoal -> 
            14, 
InterpolationPoints -> 500] // 
        Timing

The message mean that it is not possible to solve this equation on the interval [1,3). It demanded to find solution on [0,∞).

sol = 
      NISolve[y[x] - 
              3 (∫_1^(x + 1))(1/(((t + x))^2)) 
                    y[t] \[DifferentialD]t == (((E^(-1 - x)) (3 (1 + x) - 3 (E^x) (1 + 2 x) + E (1 + 3 x + 2 (x^2))))/((1 + x) (1 + 2 x))) + 
              3 (E^x) (ExpIntegralEi[-1 - 2 x] - 
                    ExpIntegralEi[-1 - x]), 
          y, {x, 1, ∞}, 
PrecisionGoal -> 
            14, 
InterpolationPoints -> 500] // 
        Timing

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

NIntegrate[Abs[sol[[2]][x] - (E^(-x))]^2, {x, 1, 40}]^(1/2)/
      NIntegrate[Abs[sol[[2]][x]]^2, {x, 1, 40}]^(1/2)

1.6

For following  integral equation the solution is (e^x)

sol = 
      NISolve[y[
                x] - (1/π) (∫_0^π) (E^(-t)) 
                    Cos[t - x Sin[ t]] y[t] \[DifferentialD]t == 
            (E^x) -  BesselJ[1, x], 
          y, {x, 0, π}, 
PrecisionGoal -> 
            15, 
PrecisionInformation -> 
            True
] // Timing

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

Sqrt[NIntegrate[Abs[(E^x) - sol[[2]][x]]^2, {x, 0, π}]]/
      Sqrt[NIntegrate[Abs[sol[[2]][x]]^2, {x, 0, π}]]

1.7

sol = 
      NISolve[z[
                x] - (∫_1^x) (t^x) BesselJ[x - 1, t] 
                  z[t] \[DifferentialD]t == 
            1 + BesselJ[x, 1] - (x^x) BesselJ[x, x], 
          z, {x, 1, 2}, 
PrecisionGoal -> 
            11, 
PrecisionInformation -> 
            True
] // Timing

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

NIntegrate[Abs[sol[[2]][x] - 1]^2, {x, 1, 2}]^(1/2)/
      NIntegrate[Abs[sol[[2]][x]]^2, {x, 1, 2}]^(1/2)

2. More complicated  Examples

2.1 ????

This solves an integral equation  with two integrals and variable limits of integration  in which first derivatives of the kernels and the solution (Sqrt[x])may  be  infinite  along the lines   y=2x , y=5x  and  x=0 correspondingly.

The solution represented in SymbolicSolution form

\(sol = 
      NISolve[w[
                x] - ∫\_x\%\(x\^2\)\@\(t - 2  x\)\ \ w[
                    t] \[DifferentialD]t + 
              4 \(∫\_0\%\(Sin[x]\)\@\(t - 5  x\)\ \ w[
                      t] \[DifferentialD]t\) == \@x + 
              1\/2\ x\ \@\(\((\(-2\) + x)\)\ x\)\ \@x\^2 - 
              1\/2\ \@\(\((\(-2\) + x)\)\ x\)\ x\^3 - 
              x\^2\ Log[\@\(-x\) + \@x] + 25\ x\^2\ Log[\@5\ \@\(-x\)] + 
              x\^2\ Log[\@\(\((\(-2\) + x)\)\ x\) + \@x\^2] - 
              25\ x\^2\ Log[\@Sin[x] + \@\(\(-5\)\ x + Sin[x]\)] - 
              5\ x\ \@Sin[x]\ \@\(\(-5\)\ x + Sin[x]\) + 
              2\ Sin[x]\^\(3/2\)\ \@\(\(-5\)\ x + Sin[x]\), 
          w, {x, 0, 1}, 
\ \ \ \ \ \ PrecisionInformation -> 
            True, 
PrecisionGoal -> 
            5, 
\ \ \ \ \ \ SolutionNorm -> 
            1, 
Error -> 
            Absolute, 
\ \ \ \ \ \ Weight -> 
            x\^\(\(-2\)/5\), 
SymbolicSolution -> 
            True] // Timing\)

This is the graph of absolute error  with weight  (x^(-2/5))

Plot[
      Evaluate[(x^(-2/5)) Abs[sol[[2]][x] - (Sqrt[x])]], {x, 0, 
        1}, PlotRange -> All, PlotStyle -> Hue[1]]

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

NIntegrate[
      Evaluate[(x^(-2/5)) Abs[sol[[2]][x] - (Sqrt[x])]], {x, 0, 
        1}]

NIntegrate[
      Evaluate[(x^(-2/5)) Abs[sol[[2]][x]]], {x, 0, 
        1}]

2.2 ????

Consider the following integral equation with three integrals

y (x) - (∫_0^(π/2))y (t) 
            Cos (x + 2 t) \[DifferentialD]t + 
        3 (∫_π^(3/2 π))y (t) (E^(-x t)) \[DifferentialD]t + 
        3 (∫_(2 π)^(3 π))y (t) (x - 
                  3 t) \[DifferentialD]t = ((9 π)/70) + ((210 (E^(-((3 π x)/2))) (1 + (E^(((π x)/2)))))/(4900 + (x^2))) + sin (70 x),  
    xϵ[0, 3  π]

with high-oscillating exact solution y(x)=sin(70x). We are seeking the solution with absolute uniform error equal to (10^(-6)).

sol = 
      NISolve[y[
                x] - (∫_0^(π/2))y[t] 
                  Cos[x + 2 t] \[DifferentialD]t + 
              3 (∫_π^(3/2 π))y[
                      t] (E^(-x t)) \[DifferentialD]t + 
              3 (∫_(2 π)^(3 π))y[
                      t] (x - 
                        3 t) \[DifferentialD]t == ((9 π)/70) + ((210 (E^(-((3 π x)/2))) (1 + (E^(((π x)/2)))))/(4900 + (x^2))) + Sin[70 x], 
          y, {x, 0, 3 π}, 
PrecisionGoal -> 
            6, 
PrecisionInformation -> True] // 
        Timing

The default value of MaxGridPoints=200 doesn't sufficient for this precision.

Actual error is rather big due to high-oscillations of the solution

(((NIntegrate[Abs[sol[[2]][x] - Sin[70 x]]^2, {x, 0, 3  π}]/NIntegrate[Abs[Sin[70 x]]^2, {x, 0, 3  π}]))^(1/2))

For achieving the required precision we must increase the value of MaxGridPoints.

sol = 
      NISolve[y[
                x] - (∫_0^(π/2))y[t] 
                  Cos[x + 2 t] \[DifferentialD]t + 
              3 (∫_π^(3/2 π))y[
                      t] (E^(-x t)) \[DifferentialD]t + 
              3 (∫_(2 π)^(3 π))y[
                      t] (x - 
                        3 t) \[DifferentialD]t == ((9 π)/70) + ((210 (E^(-((3 π x)/2))) (1 + (E^(((π x)/2)))))/(4900 + (x^2))) + Sin[70 x], 
          y, {x, 0, 3 π}, 
PrecisionGoal -> 
            6, 
PrecisionInformation -> 
            True, 
MaxGridPoints -> 500] // 
        Timing

(((NIntegrate[Abs[sol[[2]][x] - Sin[70 x]]^2, {x, 0, 3  π}]/NIntegrate[Abs[Sin[70 x]]^2, {x, 0, 3  π}]))^(1/2))

2.3

Consider  the  equation

y (x) - (∫_0^1) ker (x, t) 
              y (t) \[DifferentialD]t = (e^x), xϵ[0, 1],
ker (x, t)={({{((sh (x) sh (t - 1))/(sh (1))), 
                      0 ≤  x ≤  t}, {((sh (t) sh (x - 1))/(sh (1))), 
                      t ≤  x ≤ 1}})

with exact solution y(x)=((sh ((Sqrt[2]) (1 - x)) + 
              e  sh ((Sqrt[2 ]) x))/(sh ((Sqrt[2])))). The kernel of the equation has discontinuous partial differentials along the line x=t
and can be represented in the form.

ker[x_, t_] = 
        Which[0 ≤  x ≤  
            t, ((Sinh[x] Sinh[t - 1])/Sinh)[1], 
          t ≤  x ≤ 
            1, ((Sinh[t] Sinh[x - 1])/Sinh)[1]];

We are trying to construct  solution of the equation in the form (1) providing  relative (Subscript[L, 2])-error equal to (10^(-7)). As the kernel of the equation is non-smooth we increase the value of MaxGridPoints to 1000.

sol = 
      NISolve[y[
                x] - (∫_0^1) ker[x, t] 
                  y[t] \[DifferentialD]t == (E^x), 
          y, {x, 0, 1}, 
PrecisionGoal -> 
            7, 
MaxGridPoints -> 
            1000, 
SymbolicSolution -> True] // 
        Timing

Here we calculate the relative (Subscript[L, 2]) - error

Sqrt[
      NIntegrate[
          Evaluate[ 
            Abs[sol[[2]][x] - ((Sinh[(Sqrt[2]) (1 - x)] + E Sinh[(Sqrt[2 ]) x])/Sinh)[(Sqrt[2])]]^2], {x, 0, 1}]/
        NIntegrate[
          Evaluate[ 
            Abs[((Sinh[(Sqrt[2]) (1 - x)] + E Sinh[(Sqrt[2 ]) x])/Sinh)[(Sqrt[2])]]^2], {x, 0, 1}]]

The same equation can be solved in another way using the form (2) which much more convenient

sol = 
      NISolve[y[
                x] - (∫_x^1)(( Sinh[x] Sinh[t - 1])/Sinh)[1] 
                  y[t] \[DifferentialD]t - (∫_0^x)(( Sinh[t] 
                        Sinh[x - 1])/Sinh)[1] 
                  y[t] \[DifferentialD]t == (E^x), 
          y, {x, 0, 1}, 
SymbolicSolution -> 
            True, 
PrecisionGoal -> 14] // 
        Timing

Here we calculate the relative (Subscript[L, 2])-error

Sqrt[
      NIntegrate[
          Evaluate[ 
            Abs[sol[[2]][x] - ((Sinh[(Sqrt[2]) (1 - x)] + E Sinh[(Sqrt[2 ]) x])/Sinh)[(Sqrt[2])]]^2], {x, 0, 1}]/
        NIntegrate[
          Evaluate[ 
            Abs[((Sinh[(Sqrt[2]) (1 - x)] + E Sinh[(Sqrt[2 ]) x])/Sinh)[(Sqrt[2])]]^2], {x, 0, 1}]]

3. Two examples connected with  fractional integrals

3.1

Consider the equation

y (x) + (∫_0^x) 
            y (t) (((x - t))^(2/3)) \[DifferentialD]t == 
      (x^2) + ((27 (x^(11/3)))/220),  xϵ[0, 1]

with exact solution y (x) = (x^2).

The  default number of MaxGridPoints is insufficient for the required precision

sol = 
      NISolve[ y[
                x] + (∫_0^x) y[
                    t] (((x - t))^(2/3)) \[DifferentialD]t == 
            (x^2) + ((27 (x^(11/3)))/220), 
          y, {x, 0, 1}, 
SymbolicSolution -> 
            True, 
Error -> 
            Absolute, 
PrecisionGoal -> 6] // 
        Timing

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

√NIntegrate[
        Evaluate[Abs[sol[[2]][x] - (x^2)]^2], {x, 0, 1}]

We increase the value of MaxGridPoints to 300

sol = 
      NISolve[ y[
                x] + (∫_0^x) y[
                    t] (((x - t))^(2/3)) \[DifferentialD]t == 
            (x^2) + ((27 (x^(11/3)))/220), 
          y, {x, 0, 1}, 
SymbolicSolution -> 
            True, 
Error -> 
            Absolute, 
PrecisionGoal -> 
            6, 
MaxGridPoints -> 300] // 
        Timing

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

√NIntegrate[
        Evaluate[Abs[sol[[2]][x] - (x^2)]^2], {x, 0, 1}]

3.2

Consider more complicated example

y (x) + (∫_0^x)(((x - t))^(1/3)) 
            y (t) \[DifferentialD]t == 
      4/3 - ((2  (x^(2/3)))/(Γ (4/3) Γ (5/3))) + (x^(4/3)) - x,  xϵ[0, 1]

with exact solution  y[x]=4/3 - (2/(Gamma[4/3] Gamma[5/3])) (x^(2/3)).

sol = 
      NISolve[ y[
                x] + (∫_0^x)(((x - t))^(1/3)) 
                  y[t] \[DifferentialD]t == 
            4/3 - ((2  (x^(2/3)))/(Gamma[4/3] Gamma[5/3])) + (x^(4/3)) - 
              (x^2), y, {x, 0, 1}, 
SymbolicSolution -> 
            True, 
SolutionNorm -> ∞, 
PrecisionGoal -> 4
] // 
        Timing

Plot[
      Evaluate[Abs[sol[[2]][
              x] - (4/
                  3 - ((2  (x^(2/3)))/(Gamma[4/3] Gamma[5/3])))]], {x, 
        0, 1}, PlotRange -> All]

4. Some integral equations derived from 
functional-differential  or   difference-differential  equations

4.1

Consider the difference-differential boundary  problem

y'' (x) = q (x) y (x - 1) + f'' (x),  
    xϵ[-∞, 2]
y (-∞) = 0
y' (-∞) = 0

which is equivalent to the following integral equation

y (x) = (∫_(-∞)^(x - 1))(x - t - 1) 
            q (t + 1) y[t] \[DifferentialD]t + f (x)

We solve this problem for q (x) = (e^(x/2)) and  f (x)=(e^x)-(4/9)(e^(3/2 x - 1)). The exact solution is y (x) = (e^x).

sol = 
      NISolve[y[
              x] == (∫_(-∞)^(x - 1))(x - t - 
                      1) (E^((t + 1)/2)) 
                  y[t] \[DifferentialD]t +  
              (E^x) - (4/9) (E^(3/2 x - 1)), y, {x, -∞, 2}, 
          PrecisionGoal -> 
            11, 
InterpolationPoints -> 700] // 
        Timing

Plot[Evaluate[Abs[sol[[2]][x] - (E^x)]], {x, -11, 2}, 
      PlotRange -> All]

4.2

Consider the initial problem with retarded argument differential equation

y'' (x) = q (x) y (x/2) + f'' (x),  
    xϵ[0, 2]
y (0) = f (0)
y' (0) = f' (0)

which is equivalent to the following integral equation

y (x) = (∫_0^(x/2))(x/2 - t) q (2  t) 
            y[t] \[DifferentialD]t + f (x)

We solve this problem for  q (x)=cos(x/2) and f(x)=-(x/8) + ((9 sin (x))/8) .  
                  The exact solution is y (x) = 
            sin (x).

sol = 
      NISolve[y[
              x] == (∫_0^(x/2))(x/2 - t) Cos[t] 
                  y[t] \[DifferentialD]t  - (x/8) + ((9 Sin[x])/8), 
          y, {x, 0, 2}, 
PrecisionGoal -> 
            14, 
SolutionNorm -> ∞, 
Error -> Absolute] // Timing

Here is the graph of actual absolute error

Plot[Evaluate[Abs[sol[[2]][x] - Sin[x]]], {x, 0, 2}, 
      PlotRange -> All]

4.3

Consider a more complicated  problem with shifted and retarded arguments

y'' (x)=(1/4) y ((x/2))-(1/48) y(2 - (x/4)) p'(2 - (x/4))-(1/48) p(2 - (x/4)) y'(2 - (x/4))+f''(x),
y (0) = f (0),
y' (0)=(1/12) p (2) y (2)+f'(0).

which is equivalent to the following integral equation

y (x) = 
      2 (∫_0^(x/2))(x/2 - t) 
              y (t) \[DifferentialD]t + 
        (1/3) (∫_(2 - (x/4))^2) y (t) 
              p (t) \[DifferentialD]t + f (x).

We solve this equation for p (x) = (E^(x/2)) and f (x) = 
          1 - ((2 e)/3) + 
            (2/3) (e^(1 - (x/8))) - ((x^2)/4). The exact solution is y (x) = 1;

sol = 
      NISolve[y[x] == 
            2 (∫_0^(x/2))(x/2 - t) 
                    y[t] \[DifferentialD]t + 
              (1/3) (∫_(2 - (x/4))^2) y[
                      t] (E^(t/2)) \[DifferentialD]t + 
              1 - ((2 E)/3) + 
              (2/3) (E^(1 - (x/8))) - ((x^2)/4), 
          y, {x, -5, 4}, 
PrecisionGoal -> 14] // 
        Timing

Plot[Evaluate[Abs[sol[[2]][x] - 1]], {x, -5, -4}, 
      PlotRange -> All]

New Examples

3 Oscillatory kernel

sol = 
      NISolve[(E^(2  z)) 
                y[z] - (∫_0^1) Sin[z Log[x]]/Log[x] 
                  y[x] \[DifferentialD]x == (E^(2  z)) z -  
              ArcTan[z/2], 
          y, {z, 0, 
            1}, 
SolutionNorm -> ∞, 
PrecisionGoal -> 
            10, 
PrecisionInformation -> 
            True
] // Timing

Graph of the actual error

Plot[Abs[sol[[2]][x] - x], {x, 0, 1}, 
      PlotRange -> All]

sol = 
      NISolve[(E^(2  z)) 
                y[z] - (∫_0^1) Sin[z Log[x]]/Log[x] 
                  y[x] \[DifferentialD]x == (E^(2  z)) z -  
              ArcTan[z/2], 
          y, {z, 0, 
            1}, 
SolutionNorm -> ∞, 
PrecisionGoal -> 
            10, 
PrecisionInformation -> 
            True, 
MaxGridPoints -> 
            400, 
GaussPoints -> 
            15
] // Timing

Graph of the actual error

Plot[Abs[sol[[2]][x] - x], {x, 0, 1}, 
      PlotRange -> All]

4 Using elliptic integral (see for other examples)

(∫_0^z)(Sqrt[(( (a^2) - (x^2))/((b^2) - 
                  (x^2)))]) \[DifferentialD]x

(∫_0^(z/2))(Sqrt[(((((z + 1))^2) - (x^2))/((3^2) - 
                  (x^2)))]) \[DifferentialD]x

Solution y[x] = (1/(Sqrt[9 - (x^2)]))

sol = 
      NISolve[(Sqrt[9 - (z^2)]) 
                y[z] - (∫_0^(z/2))(Sqrt[(((z + 1))^2) - (x^2)]) 
                  y[x] \[DifferentialD]x == 
            1 - ((2 (Sqrt[9 - ((z^2)/4)]) (Sqrt[((4 + 8 z + 3 (z^2))/(36 - (z^2)))]) EllipticE[ArcSin[(z/6)], (9/(((1 + z))^2))])/(Sqrt[((4 + 8 z + 3 (z^2))/(((1 + z))^2))])), 
          y, {z, 0, 1}, 
PrecisionGoal -> 
            10, 
PrecisionInformation -> 
            True, 
SolutionNorm -> 
            4
] // Timing

Graph of the actual error

Plot[Abs[sol[[2]][x] - (1/(Sqrt[9 - (x^2)]))], {x, 0, 1}, 
      PlotRange -> All]