Autovalores molas acopladas (pontas abertas)

#include <stdio.h>
#include <math.h>

#define N 100

#define A {{1,4,8,4},{4,2,3,7},{8,3,6,9},{4,7,9,2}}


multmatrix(int n, double x[n][n],double y[n][n],double result[n][n]){
  int i,j,k;
  double soma;

  for(i=0;i<n;i++){
    for(j=0;j<n;j++){
      soma=0.0;
      for(k=0;k<n;k++)
        soma += x[i][k]*y[k][j];
      result[i][j]=soma;
    }
  }

}

imprime(int n, double x[n][n]){
  int i,j;

  for(i=0; i<n; i++){
    for(j=0; j<n; j++)
      printf("%5.2f\t",x[i][j]);
    printf("\n");
  }
  printf("\n\n");
}


QRdecomp(int n, double a[n][n], double q[n][n], double r[n][n]){
  int i,j,k;
  double u[n],prod,norm;

// ui = ai - sum_j<i (qj·ai) qj
  for(i=0;i<n;i++){

    for(k=0;k<n;k++)
      u[k] = a[k][i];

    for(j=0;j<i;j++){
//produto escalar (qj·ai)
      prod=0.0;
      for(k=0;k<n;k++)
        prod += q[k][j]*a[k][i];

      r[j][i] = prod;

      for(k=0;k<n;k++)
        u[k] -= prod*q[k][j];
    }

//qi = ui/|ui|
    norm=0.0;
    for(k=0;k<n;k++)
      norm += u[k]*u[k];
    norm = sqrt(norm);

    r[i][i] = norm;

    for(k=0;k<n;k++)
      q[k][i]=u[k]/norm;
  }

}


main(){
  int i,j,k;
  double a[N][N]={{0}};
  double q[N][N],r[N][N]={0};
  double v[N][N]={0},temp[N][N];
  double epsilon=1.0e-2;

  for(k=0;k<N;k++)
    v[k][k]=1.0;

  double kappa=6,w=2,alpha=2*kappa;

  for(i=0;i<N-1;i++){
    a[i][i]=alpha;
    a[i][i+1]= a[i][i-1] = -kappa;
  }

  i=0;   a[i][i]=alpha-kappa; a[i][i+1]=-kappa;
  i=N-1; a[i][i]=alpha-kappa; a[i][i-1]=-kappa;

  imprime(N,a);

  do{
    QRdecomp(N,a,q,r);

    multmatrix(N,r,q,a);
    multmatrix(N,v,q,r);

    for(i=0;i<N;i++){
      for(k=0;k<N;k++){
        v[i][k]=r[i][k];
        r[i][k]=0.0;
      }
    }

    // imprime(N,a);
    // imprime(N,v);

    // getchar();

  }while(fabs(a[0][1])>epsilon);

  for (int i = 0; i < N; ++i)
  {
    printf("%d %f\t", i, fabs(a[i][i]));
    for (int j = 0; j < N; ++j)
    {
      printf("%f ", v[i][j]);
    }
    printf("\n");
  }

}