Проблема с моделью изинга, находятся неправильные значения для графика

[DaniilNakhodka1995]

Если что, я создал эту тему отдельно от прошлого своего вопроса, так как прошлый решен, а теперь возник новый. Помогите, пожалуйста. Теперь ошибка с графиками, много всего с кодом перепробовал делать, на фото ниже самое лучшее, что получилось, но все равно не то. Посмотрите, пожалуйста, что может быть не так. Надо по выведенным в файлы значения построить графики в gnumplot и у меня получается так:
https://ibb.co/NNPLK7V
https://ibb.co/YZ746yn
https://ibb.co/dgT2DXH
https://ibb.co/hdJtDwP
https://ibb.co/8ryzMh6

А должно быть так:
https://ibb.co/S63mFM1
https://ibb.co/sjqyKnY
https://ibb.co/RhP3GZ1

Также, если кто знает, напишите, пожалуйста, как одновременно в gnuplot выводить 2 графика(lines и yerrorbars). На первых фото показан нужный мне вариант, а на моих все раздельно. Прописывал в gnuplot: plot "Путь к текстовому файлу" using 1:2:3 with lines/yerrorbars, но выводится только что-то одно, либо lines либо yerrorbars.

Сейчас код, с учетом всех последних изменений выглядит так:

#include <iostream>
#include <cstdlib>
#include <math.h>
#include <fstream>
#include <cmath>
#include <ctime>
using namespace std;

void array(int n, int **arr)
{
	for(int i=0; i<n; i++)
	{
		for(int j=0; j<n; j++)
		{
			arr[i][j]=1;
			cout<<arr[i][j]<<"\t";
		}
		cout<<endl;
	}
}

void sosedi(int i, int j, int n, int *s, int **arr)
{
	if(i==0)
    {
    	s[0]=arr[n-1][j];
    }
    else
    {
	    s[0]=arr[i-1][j];
    }
    if(i==n-1)
    {
    	s[1]=arr[0][j];
    }
    else
    {
	    s[1]=arr[i+1][j];
    }
    if(j==0)
    {
	    s[2]=arr[i][n-1];
    }
    else
    {
	    s[2]=arr[i][j-1];
    }
    if(j==n-1)
    {
	    s[3]=arr[i][0];
    }
    else
    {
	    s[3]=arr[i][j+1];
    }
}

double e1sp(int i, int j, int n, int **arr)
{
    double sum=0.0;
    int s[4];
    sosedi(i, j, n, s, arr);
    for(int a=0; a<4; a++)
    {
    	sum+=arr[i][j]*s[a];
    }
    return (-1)*sum;
}

double poles(int n, int **arr)
{
	double sum=0.0;
	for(int i=0; i<n; i++)
	{
		for(int j=0; j<n; j++)
		{
			sum+=e1sp(i, j, n, arr);
		}
	}
	return sum/(n*n);
}

double sreden(int n, int **arr)
{
	return poles(n, arr)/(n*n);
}

double namag(int n, int **arr)
{
	double sum=0.0;
	for(int i=0; i<n; i++)
	{
		for(int j=0; j<n; j++)
		{
			sum+=arr[i][j];
		}
	}
	return fabs(sum/(n*n));
}

double tepl(double Esum1, double Esum2, double T){
	double Esr=0.0, Esrkv=0.0;
    Esr=Esum1/10;
    Esrkv=Esum2/10;
    double t=(Esrkv-(Esr*Esr))/(T*T);
    return t;
}

double er(double *xi, double x){
    x/=10.0;
    double sum=0.0;
    for (int i=0; i<10; i++)
	{
        sum+=(xi[i]-x)*(xi[i]-x);
    }
    double m=sum*0.11111111111;
    return sqrt(m);
}

int main()
{
	int n=0;
	ofstream Ef, Mf, Cf;
	cout<<"Vvedite n:";
	cin>>n;
	cout<<endl;
	int **arr=new int *[n];
	for (int i=0; i<n; i++)
	{
		arr[i]=new int[n];
	}
	array(n, arr);
	cout<<endl;
	int i=rand()%n;
	int j=rand()%n;
	cout<<"Energia vibrannogo spina="<<e1sp(i, j, n, arr)<<endl;
	cout<<"Energia vsey sistemi="<<poles(n, arr)<<endl;
	cout<<"Namagnichennost vsey sistemi="<<namag(n, arr)<<endl;
	sreden(n, arr);
	double T, C=0.0, Esum1=0.0, Esum2=0.0;
	double Est[10], Mst[10];
	double Msum=0.0, Csum=0.0;
    Ef.open("Efull.txt");
	Mf.open("M.txt");
	Cf.open("C.txt");
	for(T = 0.01; T<4; T+=0.1)
	{
		for(int i=0; i<10; i++)
		{
			for(long int MK=0; MK<100000; MK++)
			{
				double E1=0.0, E2=0.0, p=0.0, slch=0.0;
	        	int i=rand()%n;
	            int j=rand()%n;
	            E1=e1sp(i, j, n, arr);
	            arr[i][j]=(-1)*arr[i][j];
	            E2=e1sp(i, j, n, arr);
	            if(E2>=E1)
	            {
	            	p=exp(((-1)*(E2-E1))/T);
		            slch=(double)rand()/(double)RAND_MAX;
	            	if(slch>p)
	            	{
	            		arr[i][j]=(-1)*arr[i][j];
		    		}
		    	}			 
		    }
            Est[i]=poles(n, arr);
            Esum1+=Est[i];
            Esum2+=(Est[i]*Est[i]);
            Mst[i]=namag(n, arr);
            Msum+=Mst[i];
		}
		cout<<"T="<<T<<endl;
		cout<<"E(T)="<<poles(n, arr)<<endl;
		Ef<<T<<"\t"<<poles(n, arr)<<"\t"<<er(Est, Esum1)<<endl;
		cout<<"M(T)="<<namag(n, arr)<<endl;
		Mf<<T<<"\t"<<namag(n, arr)<<"\t"<<er(Mst, Msum)<<endl;
		C=tepl(Esum1, Esum2, T);
		cout<<"C(T)="<<C<<endl;
		Cf<<T<<"\t"<<C<<endl;
	}
    for(int i=0; i<n; i++)
	{
		delete[] arr[i];	
	} 
    delete[] arr;
    Ef.close();
    Mf.close();
    Cf.close();
}

Массив специально генерируется только из единиц со знаком +, без -1, так советовал препод. Однако что с -1, что без -1, разницы особой нет, проверял.

[BDA]

Вот так выводит линии с барами, но ошибка очень большая (может быть, неверный N ввел):

plot "Efull.txt" using 1:2 with lines linetype 1, "" using 1:2:3 with errorbars linetype 1