function driver_bcg, N 
  ;N=(dindgen(1)+1.00)*4096.00
  ;make an array holding all locations of "-1" right off the main diagonal
  ;this will go through and remove locations which are multiples of N so that 
  ;the value 0 can take these places
  offDiagi=dindgen(N^2-1.0)+1.0D0
  offDiagMinus=((dindgen(N-1.0)+1.0D0)*N)-1.0D0
  offDiagi(offDiagMinus) = -1.0D0
  bad = where(offDiagi eq -1.0D0, COMPLEMENT = good)
  offDiagi = offDiagi(good)
  offDiagj=dindgen(N^2-1.0)
  offDiagj = offDiagj(good)
  offDiag = make_array(1,N^2-N, /DOUBLE, VALUE=-1);

  ;4 will be across the diagonal, or where i=j
  si=dindgen(N^2)
  sj=dindgen(N^2)
  s = make_array(1,N^2, /DOUBLE, VALUE=4)

  ;get the furthest off-diagonal locations, with all -1
  negOnei=dindgen(N^2-N)+N
  negOnej=dindgen(N^2-N)

  ;make the i and j location arrays, using the locations of the values "4", 
  ;and "-1"
  i=[si,negOnei,negOnej, offDiagi, offDiagj]
  j=[sj,negOnej,negOnei, offDiagj, offDiagi]
  si=0
  negOnei=0
  negOnej=0
  offDiagi=0
  offDiagj=0  

  ;make the array containing all the values  
  vals=[[s],[offDiag],[offDiag], [offDiag], [offDiag]]

  s=0
  offDiag=0

  ;create a sparse matrix using the ith and jth locations, and the values
  ;at those locations
  A=sprsin(i,j,vals,N^2) ;create sparse matrix
  i=0
  j=0
  vals=0

  ;set up B
  h = 1.0D0/(N+1.0D0)
  x = (dindgen(N)+1.0D0)*h
  y = x

  ;create the matrix B
  ;use CMReplicate, created by Craig B.Markwardt, to replicate an array 
  ;to create a full grid like ndgrid in matlab
  X2 = cmreplicate(x, [N,1])
  X2 = transpose(X2)
  Y2 = cmreplicate(y, [N,1]) 
  F = (-2.0D0*!PI^2.0D0)*(cos(2.0D0*!PI*X2)*(sin(!PI*Y2))^2.0D0 + (sin(!pi*X2))^2.0D0*cos(2.0D0*!pi*Y2))
  F = reform(F, N^2, 1)
  b = h^2.0 * F
  F=0
  
  start = SYSTIME(1) ;start counting how long it takes
  init=make_array(N^2, 1, /double, value=0) ;initialize an array
 ;conjugate gradient method
  u1 = linbcg(A, b, init, /DOUBLE, ITOL=1, TOL=1.0e-6, ITER=itNum, ITMAX=99999) ;call biconjugate
      ;gradient method
  Uint = reform(u1, N, N)
  stop = SYSTIME(1) ;stop counting time

  print, 'N:', N
  print, 'Time for Iterative Biconjugate Gradient:', stop-start
  print, 'Number of Iterations:', itNum
  u1=0
  
  ;create U
  x=[0, x, 1]
  y=[0, y, 1]
  X2=cmreplicate(x, [N+2,1])
  X2=transpose(X2)
  Y2=cmreplicate(y, [N+2,1])
  U=make_array(N+2,N+2,/double, value=0.0D0) 
  U[1, 1] =Uint
  
  ;calculate the true value
  Utrue= (sin(!pi*X2))^2.0D0 * (sin(!pi*Y2))^2.0D0
  X2=0
  Y2=0

  ;find the error
  Err = U-Utrue
  print, 'u-uh:', MAX(ABS(Err)) ;calculate the infinity norm of the Error

;GRAPH
imopen,'jpg', fn='bgm_graph', color=39
SURFACE, U,X,Y,charsize=2, background=255, color=0
imclose

imopen,'jpg', fn='err_graph', color=39
SURFACE, Err, X, Y, charsize=2, background=255, color=0
imclose
end