Verilog-A convergence issues

`include "constants.vams"
`include "disciplines.vams"

module test_func(p,n);

	inout p,n;
	electrical p,n;

	//files for reading/writing
	integer mcd, fp;

	real data, v_present, v_step;

	real preisach_density[0:256];// all elements should add up to 1
	real domain_states[0:256];// all elements should be -1 or +1 (depending on polarity)
	real polarization;//polarization (between -Pr and +Pr)
	
	real polarization_prev, polarization_next;
	real time_prev, time_next;

	integer i, j, k;//used in loops

	integer rows, columns;//describes area scanned from applied voltage
	integer rows_total, columns_total;//gives number of rows and columns in the grid ("rows" and "columns" should not exceed these values)
	
	analog begin	

		@(initial_step) begin

			//read data from csv file and transfer it to verilog-A array variable
			//data stored in csv file should be preisach distribution
			//csv file should have the same number of elements as arrays "preisach_density" and "domain_states"
			mcd=$fopen("test6.csv","r");

			while (!$feof(mcd)) begin

				$fscanf(mcd,"%f", preisach_density[i]);
				domain_states[i] = 1;

				i = i+1;
      			end

		//****************** v_step should be changed, should equal resolition of grid spacing
		v_step = 0.1; //arbitrary

		rows_total = 1.6/v_step;//the maximum anticipated applied input voltage divided by the resolution of the grid cells
		columns_total = 1.6/v_step;//the minimum anticipated applied input voltage divided by the resolution of the grid cells

		polarization_prev = 0;

		time_prev = 0;
		time_next = 0;
		
		end

		if (V(p,n) < 0) begin

			rows = V(p,n)/(-v_step);

			for(j = 0; j < columns_total; j=j+1) begin

				for(k=0; k<rows; k=k+1) begin

					domain_states[k*columns_total + j] = -1;

				end
			end
		end

		if (V(p,n) > 0) begin

			columns = V(p,n)/(v_step);

			//value compared to j should equal number of columns in array
			for (j=0; j<columns * rows_total ; j=j+1) begin

				domain_states[j] = 1;

			end
		end

		polarization = 0;
		for (i=0; i<columns_total*rows_total; i = i+1) begin

			polarization = polarization + preisach_density[i]*domain_states[i];

		end

		
		polarization_prev = polarization_next;
		polarization_next = polarization;
		time_prev = time_next;
		time_next = $abstime;

		I(p,n) <+ (polarization_next - polarization_prev)/(time_next - time_prev);

	end

endmodule