..................の間がプログラム本体です
program 10
..................
#include<stdio.h>
#define NSIZ 4
/* 行列演算の応用
   Gauss の消去法(ピボット選択なし)
   a x = b の解 x を求める */

main()
{
   double a[NSIZ][NSIZ],b[NSIZ],x[NSIZ];
   double awork[NSIZ][NSIZ],bwork[NSIZ],check[NSIZ];
   double dd,xx;
   int i,j,k;

/*  ファイルから行列aに値を読み込む */
  for( i=0; i<NSIZ; ++i){
    for( j=0; j<NSIZ; ++j){
      scanf("%lf",&a[i][j]);
    }
  }

/*  ファイルからベクトルbに値を読み込む */
  for( i=0; i<NSIZ; ++i){
      scanf("%lf",&b[i]);
  }

/*  作業用の行列aworkとベクトルbworkにaとbをコピー */
  for( i=0; i<NSIZ; ++i){
    for( j=0; j<NSIZ; ++j){
      awork[i][j]=a[i][j];
    }
  }
  for( i=0; i<NSIZ; ++i){
     bwork[i]=b[i];
  }

/*  前進消去 */
  for( i=1; i<NSIZ; ++i){
    for( j=i; j<NSIZ; ++j){
       dd=awork[j][i-1]/awork[i-1][i-1];
       awork[j][i-1]=0;
       bwork[j] -= dd*bwork[i-1];
       for( k=i; k<NSIZ; ++k){
          awork[j][k] -= dd*awork[i-1][k];
       }
     }
   }

/*  消去結果の出力(右上三角行列になったことを確認) */
  printf("result of forward elimination\n");
  printf("matrix awork =\n");
  for( i=0; i<NSIZ; ++i){
    for( j=0; j<NSIZ; ++j){
      printf("  %f",awork[i][j]);
    }
    printf("\n");
  }
  printf("vector bwork =\n");
  for( i=0; i<NSIZ; ++i){
    printf("  %f",bwork[i]);
  }
  printf("\n");

/*  後退代入 */
  for( i=NSIZ-1; i>=0; --i){
    xx=0;
    for( j=i+1; j<NSIZ; ++j){
       xx += x[j]*awork[i][j];
    }
    x[i] = (bwork[i]-xx)/awork[i][i];
  }

/*  代入結果の出力 */
  printf("result of backward substitution\n");
  printf("vector x =\n");
  for( i=0; i<NSIZ; ++i){
    printf("  %f",x[i]);
  }
  printf("\n");

/*  結果のチェック( a x を計算) */
  for( i=0; i<NSIZ; ++i){
    check[i] = 0;
    for( j=0; j<NSIZ; ++j){
       check[i] += a[i][j]*x[j];
    }
  }

/*  チェック結果の出力(a x - b = 0であることを確認) */
  printf("result of consistency check\n");
  printf("a x - b =\n");
  for( i=0; i<NSIZ; ++i){
    printf("  %f",check[i]-b[i]);
  }
  printf("\n");

}
..................
program 11
..................
#include<stdio.h>
#include<math.h>
#define NMAX 30
/* 関数サブプログラムの利用
 2分法による非線形方程式の解 */

double func(double);

main()
{
  double x1,x2,xtemp;
  int irepeat;

/* 解探索の範囲を指定 */
  scanf("%lf%lf",&x1,&x2);

/* 2分法を最大NMAX回くりかえす */
  for(irepeat=1; irepeat<=NMAX; ++irepeat){

/*  x1とx2の差が与えた精度以内なら終了 */
    if( fabs((x1-x2)/x1) <= 1e-10){
       printf("solution is between %17.12f and %17.12f\n",x1,x2);
       return;
    }

/*  func(x1)とfunc(x2)の符号が逆なら範囲を二分する */
    else if(func(x1)*func(x2) < 0){
       xtemp = (x1+x2)/2;
       if(func(x1)*func(xtemp) < 0){
         x2=xtemp;
       }
       else{
         x1=xtemp;
       }
       printf("%d step   x1= %17.12f   x2= %17.12f\n",irepeat,x1,x2);
    }

/*  func(x1)とfunc(x2)が同符号なら、この範囲に解がないと判断して終了 */
    else{
      printf("no solution is found between %17.12f and %17.12f\n",x1,x2);
      return;
    }
  }

/*  NMAX回のくりかえしで解に収束しなかったときは終了 */
  printf("not converged after %d steps\n",NMAX);
  printf("current values of x are %17.12f and %17.12f\n",x1,x2);

}

/*  解を得たい関数 */
double func(double x)
{
  return pow(x,3) + 3*x*x + 5*x + 1;
}