(*
	Gauss Siedel Method
	Roll no 4167441150
*)
GSiedel[A0_, n0_, x0_, xold0_, it0_, b0_] :=
  Module[{A = A0, n = n0, x = x0, xold = xold0, it = it0, b = b0},
   Upp = ConstantArray[0, {n, n}];
   For[i = 1, i <= n, i++,
    For[j = i + 1, j <= n, j++,
      Upp[[i, j]] = A[[i, j]];
      ];
    ];(* Upper triangular matrix *)

   Dia = ConstantArray[0, {n, n}];
   For[i = 1, i <= n, i++,
    Dia[[i, i]] = A[[i, i]];
    ]; (* Diagonal matrix *)

   Low = ConstantArray[0, {n, n}];
   For[i = 1, i <= n, i++,
    For[j = 1, j < i, j++,
      Low[[i, j]] = A[[i, j]];
      ];
    ]; (* Lower triangular matrix *)

   h = -(Inverse[Low + Dia]).Upp;
   c = (Inverse[Low + Dia]).b;

   (* xk+1=h.xk+c *)
   Print["\n\n-------------------OUTPUT-----------------\n"];

   For[i = 1, i <= it, i++,
    xnew = h.xold + c;
    Print["Solution at iteration #", i, " :\n"];
    For[j = 1, j <= n, j++,
     Print[x[[j, 1]], "= ", N[xnew[[j, 1]]]];
     ];
    Print["\n\n"];
    xold = xnew;
    ];
   ];

A = ({
    {10, -2, 1},
    {-2, 10, -2},
    {-2, -5, 10}
   });
b = ({
    {9},
    {12},
    {18}
   });
x = ({
    {x1},
    {x2},
    {x3}
   });
xold = ({
    {0},
    {0},
    {0}
   });
GSiedel[A, 3, x, xold, 5, b];

NRM[a0_, n_] :=
  Module[{x0 = N[a0]},
   it = 1;
   xold = 1;
   xnew = 2;
   Print["\n\n----------------------RESULT----------------------\n"]
    While[Abs[xold - xnew] != 0,
     x1 = x0 - f[x0]/f'[x0];
     Print["Root at iteration #", it++, " is: ", NumberForm[x1, 16],
      "\t Error= ", Abs[xold - xnew]];
     xold = N[x0];
     xnew = N[x1];
     x0 = x1;
     ];
   Print["\n\n Hence the root of the given equation is ",
    NumberForm[x1, 16]];
   ];
f[x_] := -0.9 x^2 + 1.7 x + 2.5;
NRM[0.5, 10];

(*
Euler's Method
Roll no 4167441150
*)
Euler[x_, y_, n_, h_] :=
  Module[{x0 = N[x], y0 = N[y]},
   Print["\n\n------------------OUTPUT--------------------\n"];
   For[i = 1, i <= n, i++,
    F = f[x0, y0];   (* y'=f(xi,yi) *)
    x0 = x0 + h;
    y0 = y0 + F*h;
    Print["Iteration #", i, ":\t 'x'= ", x0, "\t 'y'= ", y0];
    ];
   ];
f[x_, y_] := (1 + 4*x)*y^0.5;
Euler[0, 1, , 100, 0.01];