Как нарисовать траекторию полета снаряда с учетом сопротивления воздуха?

Требуется отрисовать траекторию полета снаряда с учетом силы сопротивления со стороны воздуха, которая пропорциональна квадрату скорости снаряда. Я составил рекуррентный цикл, в котором рассчитываю значения координат в каждый момент времени и меняю скорость в конце каждой итерации. Вроде с идейной точки зрения все верно, но возникает проблема: рисунок адекватен только при dt=0.1 и при малых (~0.001) значениях c. Если присвоить dt значение 0.01, то начинается какой-то ад, хотя по сути должна только вырасти точность отрисовки. Или если увеличить коэффициент c до 0.1, то траектория либо резко обрывается, либо вообще не начинается. Не понимаю, в чем дело. Важно: в целом идея программы должна остаться той же, т.е. массив из координат и рекуррентный расчет всех параметров.

Прилагаю код.

#include <windows.h>
#include "resource.h"
#include "math.h"

using namespace std;

int x1, y11, x2, y2;
int x4[465], y4[288];

int WINAPI DlgProc(HWND hDlg, WORD wMsg, WORD wParam, DWORD)
{
    PAINTSTRUCT ps;
    if (wMsg == WM_CLOSE || wMsg == WM_COMMAND && wParam == IDOK) {
        EndDialog(hDlg, 0);
    }
    else
        if (wMsg == WM_INITDIALOG) {

            RECT rc;
            GetClientRect(hDlg, &rc);
            int dx = rc.right - rc.left;
            int dy = rc.bottom - rc.top;
            cout << "dx=" << dx << " dy=" << dy << endl;

        }
        else
            if (wMsg == WM_PAINT) {
                BeginPaint(hDlg, &ps);

                HPEN hPen = (HPEN)CreatePen(PS_SOLID, 2, RGB(0, 255, 0));
                HPEN hOldPen = (HPEN)SelectObject(ps.hdc, hPen);

                POINT ptOld;
                //MoveToEx(ps.hdc, 1, 1, &ptOld);   


                hPen = (HPEN)CreatePen(PS_SOLID, 2, RGB(0, 0, 0));
                hOldPen = (HPEN)SelectObject(ps.hdc, hPen);
                MoveToEx(ps.hdc, x4[0], 288 - y4[0], &ptOld);

                for (int i = 1; i < 465; i++) {
                    LineTo(ps.hdc, x4[i], 288 - y4[i]);
                    if ((x4[i] > 465) || (y4[i] < 0) || (x4[i] < 0) || (y4[i] > 288)) {
                        break;
                    }
                }

                SelectObject(ps.hdc, hOldPen);
                DeleteObject(hPen);
                EndPaint(hDlg, &ps);

            }
    return 0;
}
void main()
{



    /* Ввод параметров задачи: 1   1    465    288    */
    cout << "Please, enter 4 coords:\n" << flush;
    //  cin >> x1 >> y11 >> x2 >> y2;
    //  Или так:
    x1 = 1;  y11 = 1;
    x2 = 465;  y2 = 288;
    cout << "x1 = " << x1 << "\ny11 = " << y11 << "\nx2 = " << x2 <<
        "\ny2 = " << y2 << "\n" << flush;


    double vx0, vy0;
    double dt;
    double c;
    double x0, y0;

    cout << "input the ratio between the drag coefficient and the mass of the object";
    cin >> c;

    cout << "input initial velocities: vx,vy:";
    cin >> vx0 >> vy0;

    cout << "input starting position(x0,y0):";
    cin >> x0 >> y0;
    double a = 1;
    x4[0] = x0;
    y4[0] = y0;
    dt = 0.01;
    for (int i = 1; i < 465; i++)
    {
        x4[i] = x4[i - 1] + vx0 * dt;
        y4[i] = y4[i - 1] + vy0 * dt;
        a = vy0 / vx0;
        vx0 = vx0 - dt * c * (vx0 * vx0 + vy0 * vy0) * cos(atan(a));
        vy0 = vy0 - dt * 9.81 - dt * c * (vx0 * vx0 + vy0 * vy0) * sin(atan(a));
    }




    DialogBox(NULL, MAKEINTRESOURCE(IDD_DIALOG1), NULL, (DLGPROC)DlgProc);
}