C / C++ / MFC

#include <afxwin.h>
#include <afxcmn.h>
#include <math.h>
#define M 200
#define IDC_BUTTON 501
#define maxInput 11
#define nMenuItems 8
	
extern double parse(CString,int,double [],int []);

class CCode10 : public CFrameWnd
{
protected:
	typedef struct
	{
		double x,y,p,z;			// components
	} PT;
	PT *pt;
	typedef struct
	{
		CString item,label;		// input string and its label
		CPoint hm;				// home coordinates
		CEdit ed;				// edit box
		CRect rc;				// rectangular region
	} INPUT;
	INPUT input[maxInput+1];
	typedef struct
	{
		CString item;			// menu item
		CPoint hm;				// home coordinates
		CRect rc;				// rectangular region
	} MENU;
	MENU menu[nMenuItems+1];
	typedef struct
	{
		CRect rc;				// rectangular region
		CPoint hm,end;			// starting and end coordinates
	} CURVE;
	CURVE curve;
	CListCtrl table;
	CFont Arial80;
	CButton btn;
	int m,idc,nInputItems,fMenu,fStatus;
public:
	CCode10();
	~CCode10();
	void ODE1RK2(),ODE1RK4(),ODE1System();
	void ODE1Taylor(),ODE1AB(),ODE2toODE1System();
	void ODE2FD1(),ODE2FD2(),SolveSLE(double **,double *);
	void DrawCurve(),ShowTable(),Clear(CRect);
	afx_msg void OnPaint(),OnLButtonDown(UINT,CPoint);
	afx_msg void OnButton();
	DECLARE_MESSAGE_MAP()
};

class CMyWinApp : public CWinApp
{
public:
	virtual BOOL InitInstance();
};</code>


and the c++ code:

<code>#include "Code10.h"
CMyWinApp MyApplication;
BOOL CMyWinApp::InitInstance()
{
    CCode10* pFrame = new CCode10;
    m_pMainWnd = pFrame;
    pFrame->ShowWindow(SW_SHOW);
    pFrame->UpdateWindow();
    return TRUE;
}

BEGIN_MESSAGE_MAP(CCode10,CFrameWnd)
	ON_WM_PAINT()
	ON_WM_LBUTTONDOWN()
    ON_BN_CLICKED(IDC_BUTTON,OnButton)
END_MESSAGE_MAP()

CCode10::CCode10()
{	
	Create(NULL,"Code10: Ordinary Differential Equations",
		WS_OVERLAPPEDWINDOW,CRect(0,0,800,635),NULL);
	Arial80.CreatePointFont(80,"Arial");
	pt=new PT [M+1];
	fMenu=0; fStatus=0; idc=301;
	curve.hm=CPoint(330,320); curve.end=CPoint(750,570);
	curve.rc=CRect(curve.hm.x-10,curve.hm.y-10,
		curve.end.x+30,curve.end.y+10);
	menu[1].item="ODE 1: Taylor (Order 3)"; 
	menu[2].item="ODE 1: Runge-Kutta 2";
	menu[3].item="ODE 1: Runge-Kutta 4"; 
	menu[4].item="ODE 1: Adams-Bashforth";
	menu[5].item="ODE 1: System"; 
	menu[6].item="ODE 2: Initial Value Problem";
	menu[7].item="ODE 2: Finite-Difference 1"; 
	menu[8].item="ODE 2: Finite-Difference 2";
	for (int i=1;i<=nMenuItems;i++)
	{
		menu[i].hm=CPoint(20,30+(i-1)*30);
		menu[i].rc=CRect(menu[i].hm.x,menu[i].hm.y,
			menu[i].hm.x+150,menu[i].hm.y+20);
	}
	input[0].hm=CPoint(200,10);
	input[0].rc=CRect(input[0].hm.x,input[0].hm.y,
		input[0].hm.x+560,input[0].hm.y+280);
	for (int i=1;i<=maxInput;i++)
		input[i].hm=CPoint(input[0].hm.x+10,input[0].hm.y+30+(i-1)*25);
}

CCode10::~CCode10()
{	
	delete pt;
}

void CCode10::OnPaint()
{
	
	CPaintDC dc(this);
	CString str;
	dc.SelectObject(Arial80);
	dc.SetBkColor(RGB(150,150,150));
	dc.SetTextColor(RGB(255,255,255));
	for (int i=1;i<=nMenuItems;i++)
	{
		dc.FillSolidRect(&menu[i].rc,RGB(150,150,150));
		dc.TextOut(menu[i].hm.x+5,menu[i].hm.y+5,menu[i].item);
	}
	dc.SetBkColor(RGB(255,255,255));
	dc.SetTextColor(RGB(100,100,100));
	dc.Rectangle(input[0].rc);
	dc.TextOut(input[0].hm.x+150,input[0].hm.y+5,menu[fMenu].item);
	if (fMenu>0)
	{
		for (int i=1;i<=nInputItems;i++)
			dc.TextOut(input[i].hm.x+10,input[i].hm.y,input[i].label);
		if (fStatus==2)
		{
			DrawCurve();
			fStatus=1;
		}
	}
	if (fMenu==2)
		dc.TextOut(input[6].hm.x+10,input[6].hm.y,input[6].label);
	if (fMenu==7)
		dc.TextOut(input[0].hm.x+300,input[0].hm.y+5,
		"p(x) y\"+q(x) y'+r(x) y = w(x),  x[0]<=x<=x[m]");
	if (fMenu==8)
		dc.TextOut(input[0].hm.x+300,input[0].hm.y+5,
		"p(x) y\"+q(x) y'+r(x) y = w(x),  x[0]<=x<=x[m]");
}

void CCode10::DrawCurve()
{
	CClientDC dc(this);
	
	double m1,m2,c1,c2;
	CString str;
	CPoint px,px2;
	PT top,bottom,max,min,left,right;
	Clear(curve.rc);
	left.x=pt[0].x; right.x=pt[m].x;
	left.y=pt[0].y;
	max.y=pt[0].y; min.y=pt[0].y;
	
	if (fMenu==5)
	{
		max.z=pt[0].z; 
		min.z=pt[0].z;
	}
	for (int i=0;i<=m;i++)
	{
		if (max.y<pt[i].y)
			max.y=pt[i].y;	
		if (min.y>pt[i].y)
			min.y=pt[i].y;
		if (fMenu==5)
		{
			if (max.z<pt[i].z)
				max.z=pt[i].z;
			if (min.z>pt[i].z)
				min.z=pt[i].z;
		}
	}
	top.y=max.y; bottom.y=min.y;
	if (fMenu==5)
	{
		if (max.y<max.z)
			top.y=max.z;
		if (min.y>min.z)
			bottom.y=min.z;
	}

	// Cartesian-Windows conversion coordinates
	m1=(double)(curve.end.x-curve.hm.x)/(right.x-left.x);
	c1=(double)curve.hm.x-left.x*m1;
	m2=(double)(curve.hm.y-curve.end.y)/(top.y-bottom.y);
	c2=(double)curve.end.y-bottom.y*m2;

	// Draw & label the x,y axes
	CPen pGray(PS_SOLID,1,RGB(100,100,100));
	dc.SelectObject(pGray);
	dc.SelectObject(Arial80);
	px=CPoint(m1*0+c1,m2*bottom.y+c2); dc.MoveTo(px);
	px=CPoint(m1*0+c1,m2*top.y+c2); dc.LineTo(px);
	
	px=CPoint(m1*left.x+c1,m2*0+c2); dc.MoveTo(px);
	str.Format("%.0lf",left.x); dc.TextOut(px.x,px.y,str);
	px=CPoint(m1*right.x+c1,m2*0+c2); dc.LineTo(px);
	str.Format("%.0lf",right.x); dc.TextOut(px.x-10,px.y,str);

	// draw the curve
	CPen pDark(PS_SOLID,2,RGB(50,50,50));
	dc.SelectObject(pDark);
	//str.Format("(%.1lf,%.1lf)",left.x,left.y);
	//dc.TextOut((int)(m1*left.x+c1),(int)(m2*left.y+c2),str);
	for (int i=0;i<=m;i++)
	{
		px=CPoint((int)(m1*pt[i].x+c1),(int)(m2*pt[i].y+c2));
		if (i==0)
			dc.MoveTo(px);
		else
			dc.LineTo(px);
		if (pt[i].y==top.y)
		{
			str.Format("%.3lf",top.y); 
			dc.TextOut(px.x,px.y-10,str);
		}
		if (pt[i].y==bottom.y)
		{
			str.Format("%.3lf",bottom.y); 
			dc.TextOut(px.x,px.y,str);
		}
	}
	if (fMenu==5)
	{
		CPen pLight(PS_SOLID,2,RGB(150,150,150));
		dc.SelectObject(pLight);
		for (int i=0;i<=m;i++)
		{
			px2=CPoint((int)(m1*pt[i].x+c1),(int)(m2*pt[i].z+c2));
			if (i==0)
			{
				dc.MoveTo(px2);
				dc.TextOut(px2.x+15,px2.y-10,"z=g(x)");
			}
			else
				dc.LineTo(px2);
		}
	}
}

void CCode10::OnLButtonDown(UINT nFlags,CPoint pt)
{
	int i,k;
	for (k=1;k<=nMenuItems;k++)
		if (menu[k].rc.PtInRect(pt))
		{
			fMenu=k; fStatus=1;
			switch(k)
			{
				case 1: 
					nInputItems=7; 
					input[1].label="y'=g(x,y)";
					input[2].label="x[0]"; 
					input[3].label="x[m]";
					input[4].label="h"; 
					input[5].label="y[0]"; 
					input[6].label="y'' (u=y')";
					input[7].label="y'''(v=y'')";
					break;
				case 2: 
					nInputItems=6; 
					input[1].label="y'=g(x,y)"; 
					input[2].label="x[0]"; 
					input[3].label="x[m]";
					input[4].label="h"; 
					input[5].label="y[0]"; 
					input[6].label="r (0<r<=1)";
					break;
				case 3: 
					nInputItems=5; 
					input[1].label="y'=g(x,y)"; 
					input[2].label="x[0]"; 
					input[3].label="x[m]";
					input[4].label="h"; 
					input[5].label="y[0]"; 
					break;
				case 4: 
					nInputItems=5; 
					input[1].label="y'=g(x,y)"; 
					input[2].label="x[0]"; 
					input[3].label="x[m]";
					input[4].label="h"; 
					input[5].label="y[0]"; 
					break;
				case 5: 
					nInputItems=7; 
					input[1].label="y'=f(x,y,z)"; 
					input[2].label="z'=g(x,y,z)"; 
					input[3].label="x[0]"; 
					input[4].label="y[0]";	
					input[5].label="z[0]";
					input[6].label="h"; 
					input[7].label="x[m]";
					break;
				case 6: 
					nInputItems=6; 
					input[1].label="y\"=f(x,y,z=y')"; 
					input[2].label="x[0]"; 
					input[3].label="y[0]"; 
					input[4].label="z[0]=y'[0]";
					input[5].label="h"; 
					input[6].label="x[m]";
					break;
				case 7: 
					nInputItems=9; 
					input[1].label="p(x)"; 
					input[2].label="q(x)"; 
					input[3].label="r(x)"; 
					input[4].label="w(x)";
					input[5].label="h"; 
					input[6].label="x[0]"; 
					input[7].label="x[m]";
					input[8].label="y[0]"; 
					input[9].label="y[m]";
					break;
				case 8: 
					nInputItems=9; 
					input[1].label="p(x)"; 
					input[2].label="q(x)"; 
					input[3].label="r(x)"; 
					input[4].label="w(x)";
					input[5].label="h"; 
					input[6].label="x[0]";
					input[7].label="x[m]"; 
					input[8].label="y'[0]";
					input[9].label="y'[m]";
					break;
			}
			for (i=1;i<=maxInput;i++)
				input[i].ed.DestroyWindow();
			for (i=1;i<=nInputItems;i++)
				input[i].ed.Create(WS_CHILD| WS_VISIBLE| WS_BORDER, 
					CRect(input[i].hm.x+100,input[i].hm.y,
					input[i].hm.x+520,input[i].hm.y+20),this,idc++);
			InvalidateRect(input[0].rc);
			Clear(curve.rc);
			table.DestroyWindow();
			btn.DestroyWindow();
			btn.Create("Compute",WS_CHILD| WS_VISIBLE| BS_DEFPUSHBUTTON,
				CRect(CPoint(input[0].hm.x+10,input[0].hm.y+5),
				CSize(100,20)),this,IDC_BUTTON);
		}
}

void CCode10::ODE1RK2()
{
	int i,psi[6];
	double h,r,psv[6],tmp,max;
	double k1,k2;
	pt[0].x=atof(input[2].item); 
	pt[0].y=atof(input[5].item);
	h=atof(input[4].item); 
	tmp=atof(input[3].item); 
	m=(tmp-pt[0].x)/h; m=((m<M)?m:M);
	max=pt[0].x+(double)m*h;
	tmp=(tmp<max)?tmp:max;
	pt[m].x=tmp;
	r=atof(input[6].item);
	psi[1]=23; psi[2]=24;
	for (i=0;i<=m;i++)
	{
		psv[1]=pt[i].x;
		psv[2]=pt[i].y;
		k1=h*parse(input[1].item,2,psv,psi);
		psv[1]=pt[i].x+r*h;
		psv[2]=pt[i].y+r*k1;
		k2=h*parse(input[1].item,2,psv,psi);
		if (i<m)
		{
			pt[i+1].y=pt[i].y+(1-1/(2*r))*k1+1/(2*r)*k2;
			pt[i+1].x=pt[i].x+h;
		}
	}
}

void CCode10::ODE1RK4()
{
	
	int i,psi[6];
	double h,psv[6],tmp,max;
	double k1,k2,k3,k4;
	pt[0].x=atof(input[2].item); 
	pt[0].y=atof(input[5].item);
	h=atof(input[4].item); 
	tmp=atof(input[3].item); 
	m=(tmp-pt[0].x)/h; m=((m<M)?m:M);
	max=pt[0].x+(double)m*h;
	tmp=(tmp<max)?tmp:max;
	pt[m].x=tmp;
	psi[1]=23; psi[2]=24;
	for (i=0;i<=m;i++)
	{
		psv[1]=pt[i].x;
		psv[2]=pt[i].y;
		k1=h*parse(input[1].item,2,psv,psi);

		psv[1]=pt[i].x+h/2;
		psv[2]=pt[i].y+k1/2;
		k2=h*parse(input[1].item,2,psv,psi);

		psv[1]=pt[i].x+h/2;
		psv[2]=pt[i].y+k2/2;
		k3=h*parse(input[1].item,2,psv,psi);

		psv[1]=pt[i].x+h;
		psv[2]=pt[i].y+k3;
		k4=h*parse(input[1].item,2,psv,psi);
		if (i<m)
		{
			pt[i+1].y=pt[i].y+(k1+2*k2+2*k3+k4)/6;
			pt[i+1].x=pt[i].x+h;
		}
	}
}

void CCode10::ShowTable()
{
	CString str;
	CPoint hTable=CPoint(20,310);
	CRect rcTable=CRect(hTable.x,hTable.y,hTable.x+280,hTable.y+290);
	table.DestroyWindow();
    table.Create(WS_VISIBLE | WS_CHILD | WS_DLGFRAME | LVS_REPORT
        | LVS_NOSORTHEADER,rcTable,this,idc++);
	table.InsertColumn(0,"i",LVCFMT_CENTER,25);
	table.InsertColumn(1,"x",LVCFMT_CENTER,70);
	table.InsertColumn(2,"y=f(x)",LVCFMT_CENTER,70);
	if (fMenu==4 || fMenu==5)
		table.InsertColumn(3,((fMenu==4)?"p":"z"),LVCFMT_CENTER,70);
	for (int i=0;i<=m;i++)
	{
		str.Format("%d",i); table.InsertItem(i,str,0);
		str.Format("%lf",pt[i].x); table.SetItemText(i,1,str);
		str.Format("%lf",pt[i].y); table.SetItemText(i,2,str);
		if (fMenu==4 && i>3)
		{
			str.Format("%lf",pt[i].p); 
			table.SetItemText(i,3,str);
		}
		if (fMenu==5)
		{
			str.Format("%lf",pt[i].z); 
			table.SetItemText(i,3,str);
		}
	}
}

void CCode10::OnButton()
{
	if (fStatus==1)
	{
		for (int i=1;i<=nInputItems;i++)
			input[i].ed.GetWindowText(input[i].item);
		switch(fMenu)
		{
			case 1:
				ODE1Taylor(); fStatus=2; break;
			case 2:
				ODE1RK2(); fStatus=2; break;
			case 3:
				ODE1RK4(); fStatus=2; break;
			case 4:
				ODE1AB(); fStatus=2; break;
			case 5:
				ODE1System(); fStatus=2; break;
			case 6:				
				ODE2toODE1System(); fStatus=2; break;
			case 7:				
				ODE2FD1(); fStatus=2; break;
			case 8:				
				ODE2FD2(); fStatus=2; break;
		}
	}
	if (fStatus=2)
	{
		ShowTable();			
		InvalidateRect(curve.rc);
	}
}

void CCode10::ODE1System()
{
	
	int i,psi[6];
	double h,psv[6],tmp,max;
	double k1,k2,k3,k4;
	double K1,K2,K3,K4;
	pt[0].x=atof(input[3].item); 
	pt[0].y=atof(input[4].item); 
	pt[0].z=atof(input[5].item);
	h=atof(input[6].item); 
	tmp=atof(input[7].item);
	m=(tmp-pt[0].x)/h; m=((m<M)?m:M);
	max=pt[0].x+(double)m*h;
	tmp=(tmp<max)?tmp:max;
	pt[m].x=tmp;
	psi[1]=23; psi[2]=24; psi[3]=25;
	for (i=0;i<=m;i++)
	{
		psv[1]=pt[i].x;
		psv[2]=pt[i].y;
		psv[3]=pt[i].z;
		k1=h*parse(input[1].item,3,psv,psi); 
		K1=h*parse(input[2].item,3,psv,psi);

		psv[1]=pt[i].x+h/2;
		psv[2]=pt[i].y+k1/2;
		psv[3]=pt[i].z+K1/2;
		k2=h*parse(input[1].item,3,psv,psi); 
		K2=h*parse(input[2].item,3,psv,psi);

		psv[1]=pt[i].x+h/2;
		psv[2]=pt[i].y+k2/2;
		psv[3]=pt[i].z+K2/2;
		k3=h*parse(input[1].item,3,psv,psi); 
		K3=h*parse(input[2].item,3,psv,psi);

		psv[1]=pt[i].x+h;
		psv[2]=pt[i].y+k3;
		psv[3]=pt[i].z+K3;
		k4=h*parse(input[1].item,3,psv,psi); 
		K4=h*parse(input[2].item,3,psv,psi);

		if (i<m)
		{
			pt[i+1].y=pt[i].y+(k1+2*k2+2*k3+k4)/6;
			pt[i+1].z=pt[i].z+(K1+2*K2+2*K3+K4)/6;
			pt[i+1].x=pt[i].x+h;
		}
	}
}

void CCode10::ODE1Taylor()
{
	double psv[6],tmp,max;
	int psi[6];
	int i;
	double h,u,v,w,z;

	pt[0].x=atof(input[2].item); 
	pt[0].y=atof(input[5].item);
	h=atof(input[4].item); 
	tmp=atof(input[3].item); 
	m=(tmp-pt[0].x)/h; m=((m<M)?m:M);
	max=pt[0].x+(double)m*h;
	tmp=(tmp<max)?tmp:max;
	pt[m].x=tmp;
	psi[1]=23; psi[2]=24;
	psi[3]=20; psi[4]=21; //psi[5]=22;
	for (i=0;i<=m;i++)
	{
		psv[1]=pt[i].x; 
		psv[2]=pt[i].y; 
		u=parse(input[1].item,2,psv,psi); 
		psv[3]=u; 
		v=parse(input[6].item,3,psv,psi); 
		psv[4]=v; 
		w=parse(input[7].item,4,psv,psi);
		if (i<m)
		{
			pt[i+1].y=pt[i].y+h*u+pow(h,2)/2*v+pow(h,3)/6*w;
			pt[i+1].x=pt[i].x+h;
		}
	}
}

void CCode10::ODE1AB()
{
	
	int i,psi[6];
	double h,f0,f1,f2,f3,fp,psv[6],tmp,max;
	double k1,k2,k3,k4;
	pt[0].x=atof(input[2].item); pt[0].y=atof(input[5].item);
	h=atof(input[4].item); 
	tmp=atof(input[3].item); 
	m=(tmp-pt[0].x)/h; m=((m<M)?m:M);
	max=pt[0].x+(double)m*h;
	tmp=(tmp<max)?tmp:max;
	pt[m].x=tmp;
	psi[1]=23; psi[2]=24;
	for (i=0;i<=m;i++)
	{
		if (i<3)
		{
			psv[1]=pt[i].x;
			psv[2]=pt[i].y;
			k1=h*parse(input[1].item,2,psv,psi);

			psv[1]=pt[i].x+h/2;
			psv[2]=pt[i].y+k1/2;
			k2=h*parse(input[1].item,2,psv,psi);

			psv[1]=pt[i].x+h/2;
			psv[2]=pt[i].y+k2/2;
			k3=h*parse(input[1].item,2,psv,psi);

			psv[1]=pt[i].x+h;
			psv[2]=pt[i].y+k3;
			k4=h*parse(input[1].item,2,psv,psi);

			pt[i+1].y=pt[i].y+(k1+2*k2+2*k3+k4)/6;
		}
		if (i<m)
		{
			pt[i+1].x=pt[i].x+h;
			if (i>=3)
			{
				psv[1]=pt[i-3].x;
				psv[2]=pt[i-3].y;
				f0=parse(input[1].item,2,psv,psi);

				psv[1]=pt[i-2].x;
				psv[2]=pt[i-2].y;
				f1=parse(input[1].item,2,psv,psi);

				psv[1]=pt[i-1].x;
				psv[2]=pt[i-1].y;
				f2=parse(input[1].item,2,psv,psi);

				psv[1]=pt[i].x; 
				psv[2]=pt[i].y; 
				f3=parse(input[1].item,2,psv,psi);
				pt[i+1].p=pt[i].y+h/24*(-9*f0+37*f1-59*f2+55*f3);

				psv[1]=pt[i+1].x;
				psv[2]=pt[i+1].p;
				fp=parse(input[1].item,2,psv,psi);
				pt[i+1].y=pt[i].y+h/24*(f1-5*f2+19*f3+9*fp);
			}
		}
	}
}

void CCode10::ODE2toODE1System()
{
	
	int i,psi[6];
	double psv[6];
	double h,tmp;
	double k1,k2,k3,k4;
	double K1,K2,K3,K4;
	pt[0].x=atof(input[2].item); 
	pt[0].y=atof(input[3].item); 
	pt[0].z=atof(input[4].item);
	h=atof(input[5].item); 
	tmp=atof(input[6].item); 
	m=(tmp-pt[0].x)/h; m=((m<M)?m:M);
	pt[m].x=tmp;
	psi[1]=23; psi[2]=24; psi[3]=25;
	for (i=0;i<=m;i++)
	{
		psv[1]=pt[i].x;
		psv[2]=pt[i].y;
		psv[3]=pt[i].z;
		k1=h*pt[i].z; K1=h*parse(input[1].item,3,psv,psi);

		psv[1]=pt[i].x+h/2;
		psv[2]=pt[i].y+k1/2;
		psv[3]=pt[i].z+K1/2;
		k2=h*(pt[i].z+K1/2); K2=h*parse(input[1].item,3,psv,psi);

		psv[1]=pt[i].x+h/2; 
		psv[2]=pt[i].y+k2/2;
		psv[3]=pt[i].z+K2/2;
		k3=h*(pt[i].z+K2/2); K3=h*parse(input[1].item,3,psv,psi);

		psv[1]=pt[i].x+h; 
		psv[2]=pt[i].y+k3;
		psv[3]=pt[i].z+K3;
		k4=h*(pt[i].z+K3); K4=h*parse(input[1].item,3,psv,psi);

		if (i<m)
		{
			pt[i+1].y=pt[i].y+(k1+2*k2+2*k3+k4)/6;
			pt[i+1].z=pt[i].z+(K1+2*K2+2*K3+K4)/6;
			pt[i+1].x=pt[i].x+h;
		}
	}
}

void CCode10::ODE2FD1()
{
	
	int i,j,psi[6];
	double psv[6],h,tmp;
	double **a,*b;
	double *p,*q,*r,*w;
	b=new double [M+1];
	p=new double [M+1];
	q=new double [M+1];
	r=new double [M+1];
	w=new double [M+1];
	a=new double *[M+1];
	h=atof(input[5].item);
	m=(int)(atof(input[7].item)-atof(input[6].item))/h; m=((m<M)?m:M);
	pt[0].x=atof(input[6].item); 
	pt[m].x=pt[0].x+(double)m*h;
	pt[0].y=atof(input[8].item); 
	pt[m].y=atof(input[9].item);
	psi[1]=23;
	for (i=1;i<=m-1;i++)
	{
		pt[i].x=pt[i-1].x+h;
		psv[1]=pt[i].x; 
		p[i]=parse(input[1].item,1,psv,psi);
		q[i]=parse(input[2].item,1,psv,psi);
		r[i]=parse(input[3].item,1,psv,psi);
		w[i]=parse(input[4].item,1,psv,psi);
	}
	for (i=0;i<=M;i++)
		a[i]=new double [M+1];
	for (i=1;i<=m-1;i++)
	{
		for (j=1;j<=m-1;j++)
			a[i][j]=0;
		b[i]=0;
	}
	b[1]=w[1]-(p[1]/(h*h)-q[1]/(2*h))*pt[0].y;
	b[m-1]=w[m-1]-(p[m-1]/(h*h)+q[m-1]/(2*h))*pt[m].y;
	for (i=1;i<=m-1;i++)
	{
		a[i][i]=-2*p[i]/(h*h)+r[i];
		if (i<m-1)
			a[i+1][i]=p[i+1]/(h*h)-q[i+1]/(2*h);
		if (i>1)
			a[i-1][i]=p[i-1]/(h*h)+q[i-1]/(2*h);
		if (i>1 && i<m-1)
			b[i]=w[i];
	}
	SolveSLE(a,b);
	for (i=0;i<=M;i++)
		delete a[i];
	delete p,q,r,w,a,b;
}

void CCode10::SolveSLE(double **a,double *b)
{
	int i,j,k,lo,hi;
	double m1,Sum;
	if (fMenu==7)
	{
		lo=1; hi=m-1;
	}
	if (fMenu==8)
	{
		lo=0; hi=m;
	}
	for (k=lo;k<=hi-1;k++)
		for (i=k+1;i<=hi;i++)
		{
			m1=a[i][k]/a[k][k];
			for (j=lo;j<=hi;j++)
				a[i][j] -= m1*a[k][j];
			b[i] -= m1*b[k];
		}
	for (i=hi;i>=lo;i--)
	{
    	Sum=0;
		pt[i].y=0;
		for (j=i;j<=hi;j++)
       		Sum += a[i][j]*pt[j].y;
		pt[i].y=(b[i]-Sum)/a[i][i];
	}
}

void CCode10::ODE2FD2()
{
	
	int i,j,psi[6];
	double psv[6],tmp;
	double h,alpha,beta;
	double **a,*b;
	double *p,*q,*r,*w;
	b=new double [M+1];
	p=new double [M+1];
	q=new double [M+1];
	r=new double [M+1];
	w=new double [M+1];
	a=new double *[M+1];
	for (i=0;i<=M;i++)
		a[i]=new double [M+1];
	h=atof(input[5].item);
	m=(int)(atof(input[7].item)-atof(input[6].item))/h; m=((m<M)?m:M);
	pt[0].x=atof(input[6].item); 
	pt[m].x=pt[0].x+(double)m*h;
	alpha=atof(input[8].item); beta=atof(input[9].item); 
	for (i=0;i<=m;i++)
	{
		if (i<m)
			pt[i+1].x=pt[i].x+h;
		psv[1]=pt[i].x; psi[1]=23;
		p[i]=parse(input[1].item,1,psv,psi);
		q[i]=parse(input[2].item,1,psv,psi);
		r[i]=parse(input[3].item,1,psv,psi);
		w[i]=parse(input[4].item,1,psv,psi);
	}
	for (i=0;i<=m;i++)
	{
		for (j=0;j<=m;j++)
			a[i][j]=0;
		b[i]=0;
	}
	for (i=0;i<=m;i++)
	{
		a[i][i]=-2*p[i]/(h*h)+r[i];
		if (i>0 && i<m)
			a[i][i+1]=p[i]/(h*h)+q[i]/(2*h);
		if (i<m-1)
			a[i+1][i]=p[i+1]/(h*h)-q[i+1]/(2*h);
		if (i>0 && i<m)
			b[i]=w[i];
	}
	a[0][1]=2*p[0]/m;
	a[m][m-1]=2*p[m]/(h*h);
	b[0]=w[0]+(p[0]/(h*h)-q[0]/(2*h))*2*h*alpha;
	b[m]=w[m]-(p[m]/(h*h)+q[m]/(2*h))*2*h*beta;
	SolveSLE(a,b);
	for (i=0;i<=M;i++)
		delete a[i];
	delete p,q,r,w,a,b;
}

void CCode10::Clear(CRect rc)
{
	CClientDC dc(this);
	CBrush whiteBrush(RGB(255,255,255));
	dc.FillRect(&rc,&whiteBrush);
}

/* Paint_beginner.cpp *********************************************************
Author:		Paul Watt
Date:		2/28/2002
Purpose:	Implementation file for a beginner's guide to painting with windows.
			
			Normally I like to create object oriented encapsulations of the 
			Main Window and data, however this program is very light.  There 
			are not any features that do not need to exist, in order to keep 
			the code clean, and to hopefully help the painting code stand out.

			Most of the data will be stored in global variables for clarity, not
			style.
******************************************************************************/
#include "stdafx.h"
#include "resource.h"

#define MAX_LOADSTRING 100

/* Global Variables **********************************************************/
HINSTANCE hInst;
TCHAR szTitle[MAX_LOADSTRING];
TCHAR szWindowClass[MAX_LOADSTRING];	

typedef struct tagShape
{
	RECT rect;
	UINT shapeID;
} Shape;


const UINT SHAPE_COUNT = 5;
Shape Shapes[SHAPE_COUNT]; 
UINT nCurrentIndex = 0;

bool  isRubberBand = false;
POINT ptStart;
POINT ptCurrent;

/* Forward Declarations ******************************************************/
ATOM				MyRegisterClass(HINSTANCE hInstance);
BOOL				InitInstance(HINSTANCE, int);
LRESULT CALLBACK	WndProc(HWND, UINT, WPARAM, LPARAM);
LRESULT CALLBACK	About(HWND, UINT, WPARAM, LPARAM);

/* Global *********************************************************************
Author:		Paul Watt
Date:		2/28/2002
Purpose:	Main entry point for this application.
******************************************************************************/
int APIENTRY WinMain(HINSTANCE hInstance,
                     HINSTANCE hPrevInstance,
                     LPSTR     lpCmdLine,
                     int       nCmdShow)
{
	MSG msg;
	HACCEL hAccelTable;
			//C: Initialize the global strings.
	LoadString(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
	LoadString(hInstance, IDC_PAINT_BEGINNER, szWindowClass, MAX_LOADSTRING);
			//C: Register the class for the main window of this application.
	WNDCLASSEX wcex;

	wcex.cbSize = sizeof(WNDCLASSEX); 

	wcex.style			= CS_HREDRAW | CS_VREDRAW;
	wcex.lpfnWndProc	= (WNDPROC)WndProc;
	wcex.cbClsExtra		= 0;
	wcex.cbWndExtra		= 0;
	wcex.hInstance		= hInstance;
	wcex.hIcon			= LoadIcon(hInstance, (LPCTSTR)IDI_PAINT_BEGINNER);
	wcex.hCursor		= LoadCursor(NULL, IDC_ARROW);
	wcex.hbrBackground	= (HBRUSH)(COLOR_WINDOW+1);
	wcex.lpszMenuName	= (LPCSTR)IDC_PAINT_BEGINNER;
	wcex.lpszClassName	= szWindowClass;
	wcex.hIconSm		= LoadIcon(wcex.hInstance, (LPCTSTR)IDI_SMALL);

	RegisterClassEx(&wcex);
			//C: Perform application initialization.
	if (!InitInstance (hInstance, nCmdShow)) 
	{
		return FALSE;
	}

	hAccelTable = LoadAccelerators(hInstance, (LPCTSTR)IDC_PAINT_BEGINNER);
			//C: Main message pump.
	while (GetMessage(&msg, NULL, 0, 0)) 
	{
		if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg)) 
		{
			TranslateMessage(&msg);
			DispatchMessage(&msg);
		}
	}

	return msg.wParam;
}


/* Global *********************************************************************
Author:		Paul Watt
Date:		2/28/2002
Purpose:	Initializes global variables and creates the MainWindow.
Parameters:	hInstance[in]: The HINSTANCE of this application.
			nCmdShow[in]: The display mode to set for the main window when it 
				is created.
Return:		If this function succeeds then true will be returned otherwise false.
******************************************************************************/
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
   HWND hWnd;
			//C: Store the instance handle in the global varaible.
   hInst = hInstance;
			//C: Initialize all of the shapes rectangles to -1.
   memset(&Shapes, -1, sizeof(Shapes));
			//C: Create the mainwindow.
   hWnd = CreateWindow(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, NULL, NULL, hInstance, NULL);
   if (!hWnd)
   {
			//C: The main window creation failed.
      return FALSE;
   }
			//C: Display the main window.
   ShowWindow(hWnd, nCmdShow);
			//C: Force the main window to repaint itself.
   UpdateWindow(hWnd);

   return TRUE;
}


LRESULT OnCommand (HWND hWnd, int iID, int iEvent, HWND hWndControl, bool &isHandled);
LRESULT OnLButtonDown (HWND hWnd, UINT nCtrl, UINT x, UINT y);
LRESULT OnMouseMove   (HWND hWnd, UINT nCtrl, UINT x, UINT y);
LRESULT OnLButtonUp   (HWND hWnd, UINT nCtrl, UINT x, UINT y);
LRESULT OnPaint       (HWND hWnd);
void DrawRubberBand(HWND hWnd);


/* Global *********************************************************************
Author:		Paul Watt
Date:		2/28/2002
Purpose:	Window procedure for the main window.
Parameters:	hWnd[in]: Handle to the window that the message is intended for.
			message[in]: The message being sent.
			wParam[in]: Data
			lParam[in]: Data
Return:		
******************************************************************************/
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
	switch (message) 
	{
	case WM_CREATE:
		{
			//C: Set the initial drawing mode.
			HMENU hMenu = ::GetMenu(hWnd);
			HMENU hMenuShapes = ::GetSubMenu(hMenu, 1);
			::CheckMenuRadioItem(hMenuShapes, ID_SHAPE_RECTANGLE, ID_SHAPE_CIRCLE, ID_SHAPE_RECTANGLE, MF_BYCOMMAND);
		}
	case WM_COMMAND:
		{
			int wmId    = LOWORD(wParam); 
			int wmEvent = HIWORD(wParam); 

			bool isHandled = true;
			LRESULT lResult = OnCommand(hWnd, wmId, wmEvent, (HWND)lParam, isHandled);
			if (!isHandled)
			{
				lResult = DefWindowProc(hWnd, message, wParam, lParam);
			}

			return lResult;
		}
		break;
	case WM_CHAR:
		{
			if (VK_ESCAPE == wParam && isRubberBand)
			{
			//C: Set rubberbanding to false.
				isRubberBand = false;
			//C: Erase the last rubber band rectangle that was drawn.
				DrawRubberBand(hWnd);
			}
		}
		break;
	case WM_LBUTTONDOWN:
		{
			return OnLButtonDown(hWnd, (UINT)wParam, LOWORD(lParam), HIWORD(lParam));
		}
		break;
	case WM_MOUSEMOVE:
		{
			return OnMouseMove(hWnd, (UINT)wParam, LOWORD(lParam), HIWORD(lParam));
		}
		break;
	case WM_LBUTTONUP:
		{
			return OnLButtonUp(hWnd, (UINT)wParam, LOWORD(lParam), HIWORD(lParam));
		}
		break;
	case WM_PAINT:
		{
			return OnPaint(hWnd);
		}
		break;
	case WM_DESTROY:
		{
			//C: Send a shutdown message to the message pump.
			PostQuitMessage(0);
			return 0;
		}
		break;
	}
			//C: Exit.
	return DefWindowProc(hWnd, message, wParam, lParam);
}


/* Global *********************************************************************
Author:		Paul Watt
Date:		2/28/2002
Purpose:	WM_COMMAND message handler.  Basically parses the menu commands
			and either sets the current drawing shape, or exits the window.
Parameters:	iID[in]:
			iEvent[in]:
			hWndControl[in]:
			isHandled[out]: This value can be set to indicate if this function 
				handled the message or not.
Return:		The LRESULT of this handled message will be returned.
******************************************************************************/
LRESULT OnCommand (HWND hWnd, int iID, int iEvent, HWND hWndControl, bool &isHandled)
{
			//C: Parse the menu selections.
	switch (iID)
	{			
	case ID_SHAPE_RECTANGLE:
	case ID_SHAPE_ELLIPSE:
		{
			//C: Set the drawing mode.
			HMENU hMenu = ::GetMenu(hWnd);
			HMENU hMenuShapes = ::GetSubMenu(hMenu, 1);
			::CheckMenuRadioItem(hMenuShapes, ID_SHAPE_RECTANGLE, ID_SHAPE_ELLIPSE, iID, MF_BYCOMMAND);
		}
		break;
	case IDM_EXIT: 
		{
			//C: Destroy the window in order to exit the program.
			DestroyWindow(hWnd);
		}
		break;
	default:
		{
			//C: Flag this message as unhandled.
			isHandled = false;
		}
	}

	return 0;
}


/* Global *********************************************************************
Author:		Paul Watt
Date:		2/28/2002
Purpose:	Handles the WM_LBUTTONDOWN message.  This function will start the
			rubberbanding effect when the user drags the mouse until the
			left button is released.  If the Escape button is pressed before the 
			left button is released, then the action will be cancelled. 
Parameters:	hWnd[in]: The handle of the window where the mouse button was pressed.
			nCtrl[in]: A set of flags that indicates the control keys that 
				are currently pressed.
			x[in]: The X coordinate where the left button was pressed.
			y[in]: The Y coordinate where the left button was pressed.
Return:		0 will be returned.
******************************************************************************/
LRESULT OnLButtonDown (HWND hWnd, UINT nCtrl, UINT x, UINT y)
{
			//C: Set the start point.
	ptStart.x = x;
	ptStart.y = y;
			//C: Set the rubberband mode to true.
	isRubberBand = true;
			//C: Set the current point to the start point as well.
	ptCurrent.x = x;
	ptCurrent.y = y;
			//C: Draw the initial rubber band.
	DrawRubberBand(hWnd);
			//C: Return.
	return 0;	
}


/* Global *********************************************************************
Author:		Paul Watt
Date:		2/28/2002
Purpose:	Handles the WM_MOUSEMOVE message.  This function will continue the
			rubberbanding effect when the user drags the mouse until the
			left button is released.  If the Escape button is pressed before the 
			left button is released, then the action will be cancelled. 
Parameters:	hWnd[in]: The handle of the window where the mouse button was pressed.
			nCtrl[in]: A set of flags that indicates the control keys that 
				are currently pressed.
			x[in]: The X coordinate where the left button was pressed.
			y[in]: The Y coordinate where the left button was pressed.
Return:		0 will be returned.
******************************************************************************/
LRESULT OnMouseMove   (HWND hWnd, UINT nCtrl, UINT x, UINT y)
{
			//C: If the current mode is not rubber band, then exit.
	if (!isRubberBand)
	{
		return 0;
	}
			//C: Undo the last rectangle shape that was drawn.
	DrawRubberBand(hWnd);
			//C: Update the current position.
	ptCurrent.x = x;
	ptCurrent.y = y;
			//C: Draw the next rubberband position.
	DrawRubberBand(hWnd);
			//C: Exit with success.
	return 0;		
}


/* Global *********************************************************************
Author:		Paul Watt
Date:		2/28/2002
Purpose:	Handles the WM_LBUTTONDOWN message.  This function will end the
			rubberbanding effect and add a new shape to the current set of shapes.
			If the Escape button is pressed before the left button is released, 
			then the action will be cancelled. 
Parameters:	hWnd[in]: The handle of the window where the mouse button was pressed.
			nCtrl[in]: A set of flags that indicates the control keys that 
				are currently pressed.
			x[in]: The X coordinate where the left button was pressed.
			y[in]: The Y coordinate where the left button was pressed.
Return:		0 will be returned.
******************************************************************************/
LRESULT OnLButtonUp   (HWND hWnd, UINT nCtrl, UINT x, UINT y)
{
			//C: If the current mode is not rubber band, then exit.
	if (!isRubberBand)
	{
		return 0;
	}
			//C: Undo the last rectangle shape that was drawn.
	DrawRubberBand(hWnd);
			//C: Update the dimensions of the current shape.
			//C: Get the current shape mode.
	HMENU hMenu			= ::GetMenu(hWnd);
	HMENU hShapeMenu	= ::GetSubMenu(hMenu, 1);

	if (::GetMenuState(hShapeMenu, ID_SHAPE_RECTANGLE, MF_BYCOMMAND) & MF_CHECKED)
	{
		Shapes[nCurrentIndex].shapeID = ID_SHAPE_RECTANGLE;
	}
	else
	{
		Shapes[nCurrentIndex].shapeID = ID_SHAPE_ELLIPSE;
	}
	
	Shapes[nCurrentIndex].rect.left		= ptStart.x;
	Shapes[nCurrentIndex].rect.top		= ptStart.y;
	Shapes[nCurrentIndex].rect.right	= ptCurrent.x;
	Shapes[nCurrentIndex].rect.bottom	= ptCurrent.y;
			//C: Increment the index.
	nCurrentIndex = (++nCurrentIndex)%SHAPE_COUNT;
			//C: Disable the rubberband mode.
	isRubberBand = false;
			//C: Update the window.
	InvalidateRect(hWnd, NULL, TRUE);
	UpdateWindow(hWnd);
			//C: Return with success.
	return 0;		
}


/* Global *********************************************************************
Author:		Paul Watt
Date:		2/28/2002
Purpose:	Paints the five shapes that have been set on this window.
Parameters:	hWnd[in]:
Return:		
******************************************************************************/
LRESULT OnPaint       (HWND hWnd)
{
	PAINTSTRUCT ps;
	HDC			hdc;
	hdc = ::BeginPaint(hWnd, &ps);

	UINT index;
	for (index = 0; index < SHAPE_COUNT; index++)
	{
		if (ID_SHAPE_RECTANGLE == Shapes[index].shapeID)
		{
			::Rectangle	(	
						hdc, 
						Shapes[index].rect.left, 
						Shapes[index].rect.top, 
						Shapes[index].rect.right,
						Shapes[index].rect.bottom
						);
		}
		else
		{
			::Ellipse	(	
						hdc, 
						Shapes[index].rect.left, 
						Shapes[index].rect.top, 
						Shapes[index].rect.right,
						Shapes[index].rect.bottom
						);
		}
	}

	::EndPaint(hWnd, &ps);

	return 0;	
}


/* Global *********************************************************************
Author:		Paul Watt
Date:		3/1/2002
Purpose:	Draws the rubberband shape using the dimensions of the start
			and current point global variables.
Parameters:	NONE
Return:		-
******************************************************************************/
void DrawRubberBand(HWND hWnd)
{
	HDC hdc;
			//C: Get a client DC.
	hdc = ::GetDC(hWnd);
			//C: Set the current drawing mode to XOR, this will allow us
			//   to add the rubber band, and later remove it by sending the
			//   exact same drawing command.
	::SetROP2(hdc, R2_NOT);
			//C: Select a NULL Brush into the DC so that no fill is performed.
	::SelectObject(hdc, ::GetStockObject(NULL_BRUSH));
			//C: Get the current shape mode.
	HMENU hMenu			= ::GetMenu(hWnd);
	HMENU hShapeMenu	= ::GetSubMenu(hMenu, 1);

	if (::GetMenuState(hShapeMenu, ID_SHAPE_RECTANGLE, MF_BYCOMMAND) & MF_CHECKED)
	{
		::Rectangle	(	
					hdc, 
					ptStart.x, 
					ptStart.y, 
					ptCurrent.x,
					ptCurrent.y
					);
	}
	else
	{
		::Ellipse	(	
					hdc, 
					ptStart.x, 
					ptStart.y, 
					ptCurrent.x,
					ptCurrent.y
					);
	}
			//C: Release the DC.
	::ReleaseDC(hWnd, hdc);
}