C Progam for calculating the stability regions of Mathieu's equation

Many thanks to Christian Schneider for spotting typos here!

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

int main(){
  FILE *fp;   // Prepare to print to file
  fp=fopen("mat.dat","w");
  int m,i,j;
  float e[201],d[101],alpha,beta,alpha1,mu,a,q;
  float pi=3.141592653589;
  a=0.5;
  q=0;

  for(q=-10;q<10;q+=0.02){   //Loop over the desired
    for(a=-5;a<10;a+=0.07){  //a-q region

  for(m=0;m<=248;m+=2){
  e[m]=q/((m*m*1.0)-a);}  //Set all components

  //The first seed determinants, from Maple worksheet
  d[3]=-2*e[2]*e[2]*e[0]*e[4]*e[4]*e[6]+e[2]*e[2]*e[4]*e[4]
           -2*e[4]*e[4]*e[2]*e[0]*e[6]*e[6]+2*e[2]*e[4]*e[4]*e[6]
           +e[4]*e[4]*e[6]*e[6]+2*e[2]*e[2]*e[0]*e[4]+4*e[2]*e[0]*e[6]*e[4]-
           2*e[2]*e[4]-2*e[6]*e[4]-2*e[2]*e[0]+1;

  d[2]=1-2*e[2]*e[4]-2*e[2]*e[0]+2*e[2]*e[2]*e[0]*e[4]+e[2]*e[2]*e[4]*e[4];
   
  d[1]=1-2*e[2]*e[0];
  
  d[0]=1;
  
  for(m=4; m<=100; m++){       //Here goes Strang's interation method
    alpha=e[2*m]*e[2*(m-1)];
    beta=1-alpha;
    alpha1=e[2*(m-1)]*e[2*(m-2)];
    d[m]=beta*d[m-1] - alpha*beta*d[m-2] + alpha*alpha1*alpha1*d[m-3];  }

  //Find mu, make seperate case for -a situation
  if(a>=0){  mu=acos( 1 - (d[100])*(1-cos(pi*sqrt(a)))) / (pi);}
  if(a<0){   mu=acos( 1 - (d[100])*(1-cosh(pi*sqrt(fabs(a))))) / (pi);}
  if (mu != mu){mu=0.000000;}  //If mu=nan then make it zero
 
  fprintf(fp,"%f %f  %f \n", q, a, mu);}fprintf(fp,"\n");}  
}

As a final note, we wish to add the code used for gnuplot. One run was to outline the edge of the $\mu=0$ isobar, since gnuplot does not plot contour this (starts with 0.05). The second run does the overall contouring as seen in the above figures.

For the border outline:

set data style lines
set contour base
set cntrparam levels discrete 0.001
set nosurface
set view 0,0
splot 'so2.txt'

For the inner contours:

set data style lines
set cntrparam levels 20
set contour
set nosurface
set view 0,0
splot 'so2.txt'